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ABSTRACT 


Time-Frequency  Analysis  in  Radar  Backscatter  Problems 


by 

Christopher  Joseph  McCormack 


Co-Chairs:  William  J.  Williams,  Valdis  V.  Liepa 

Time-frequency  techniques  provide  new  and  unique  insights  for  analyzing  electromag¬ 
netic  scattering  problems.  These  techniques  transform  functions  of  time  or  frequency  into 
two  dimensional  functions  of  both  time  and  frequency  to  reveal  non-stationary  charac¬ 
teristics  of  the  signal.  The  theory  developed  herein  justifies  applying  the  frequency-time 
transform  to  wide  bandwidth  signals  illuminating  stationary  targets.  The  frequency-time 
representation  of  the  return  provides  more  information  about  the  target  and  the  scattering 
than  regular  Fourier  analysis.  Along  with  the  position  of  the  scattering  centers,  frequency¬ 
time  analysis  gives  insights  on  the  target’s  composition  and  configuration.  In  addition, 
the  performance  of  these  transforms  when  applied  to  noise  are  examined  and  quantified. 
The  statistics  of  a  Cohen’s  class  time-frequency  transformation  are  derived  and  verified  nu¬ 
merically.  Applying  the  time-frequency  techniques  to  sampled  continuous-wave  radar  data 
from  a  dynamic  target  provides  insight  into  target  motion  and  generate  estimates  of  the 
target  parameters.  After  considering  a  figure  of  merit  for  evaluating  the  time-frequency 
distribution,  a  customized  kernel,  determined  using  a  genetic  algorithm,  is  used  to  improve 
the  performance  of  the  standard  spectrogram,  the  Wigner  distribution,  and  the  binomial 
distribution. 
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CHAPTER  1 


INTRODUCTION 


1.1  Analysis  Domains 

Many  different  representations  may  be  used  to  analyze  signals.  The  transformation  of  a 
signal  from  one  form  to  another  may  allow  simplification  of  the  math  (for  example,  replacing 
convolutions  with  multiplications),  or  it  may  provide  better  insight  into  the  behavior  of  a 
signal  (or  system)  by  relating  the  measured  signal  to  expected  signal  characteristics. 

Some  commonly  used  one-dimensional  transforms  include  the  Fourier  transform,  the 
Laplace  transform,  and  the  Z-transform.  Another  type  of  transform  moves  from  one  di¬ 
mension  to  multiple  dimensions.  These  multidimensional  transforms  also  include  Cohen’s 
class  of  time-frequency  distributions  [12,  32]  and  the  wavelet  transform  [6,  32]. 

1.1.1  Time  Domain 

The  time  domain  representation  for  a  signal  typically  represents  the  most  natural  form 
of  a  signal.  Conditions  of  realizability  force  any  real-world  signal  to  have  finite  starting 
and  ending  times.  This  domain  also  relates  directly  to  the  direct  measurements  that  can 
be  obtained.  The  time  domain  representation  provides  the  simplest  view  of  aperiodic  or 
non-stationary  signals. 

When  accurately  measured,  the  time  domain  of  the  signal  contains  all  the  information 
available  in  the  signal.  Other  domains  might  make  certain  relationships  more  clearly  visible, 
but  no  additional  information  is  generated. 
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1.1.2  Frequency  Domain 

Fourier  analysis  is  normally  introduced  by  starting  with  a  function  of  time,  g{t),  and 
finding  an  equivalent  function  of  frequency,  G{f).  This  uses  the  Fourier  integral  to  substi¬ 
tute  a  dependence  on  the  frequency  variable  /  instead  of  the  time  variable  t.  The  frequency 
domain  representation  emphasizes  the  periodic  aspects  of  a  signal.  For  linear  system  anal¬ 
ysis,  this  provides  a  valuable  analysis  tool,  both  for  insight  into  the  nature  of  a  system  or 
signal  and  for  evaluation  of  a  system’s  response. 

It  is  important  to  remember  that  the  frequency  domain  representation  is  a  mathematical 
concept.  Theoretically,  a  finite  time  signal  will  have  a  frequency  domain  representation  that 
is  not  band  limited.  Infinite  frequency  limits  are  not  a  problem  when  the  Fourier  domain 
is  treated  mathematically. 

1.1.3  Time-Frequency  Domain 

While  both  the  time  and  frequency  domain  representations  of  a  signal  are  complete, 
certain  relationships  are  more  apparent  in  the  time  domain  (signal  transitions,  starts,  and 
stops),  while  other  effects  are  easier  to  identify  in  the  frequency  domain  (noise  level,  signal 
bandwidth).  Time  frequency  analysis  attempts  to  maintain  both  time  and  frequency  as 
joint  variables  to  show  relationships  based  on  their  interactions. 

1.2  Scattering  Problems 

Time-fi:equency  representations  have  been  used  extensively  for  analyzing  non-stationary 
systems  where  the  signal  spectra  varies  as  a  fimction  of  time.  In  some  cases,  the  objective 
was  to  characterize  individual  components  of  the  signal.  In  other  cases,  the  overall  patterns 
in  the  time-frequency  representation  served  as  the  metric  for  making  decisions  about  the 
signal  [1,  12]. 

Most  of  the  work  using  time-frequency  techniques  has  focused  on  biological  and  biomed¬ 
ical  signals  [1,  43].  Very  little  had  been  done  to  address  the  use  of  these  transforms  for 
problems  involving  electromagnetic  scattering.  The  remainder  of  this  dissertation  explains 
the  application  of  different  time-firequency  procedures  to  several  different  types  of  scattering 
problems. 
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1.2.1  Dispersive  Phenomena 

Traditional  signal  analysis  for  signals  resulting  from  electromagnetic  scattering  assumes 
stationary,  non-dispersive  medium.  Standard  frequency  domain  techniques  cannot  properly 
identify  effects  which  vary  as  a  function  of  frequency.  This  is  important  for  most  non- 
idealized  scattering  targets.  The  dispersive  natiire  of  non-metallic  materials  and  complex 
scattering  modes  tends  to  smear  the  transformed  response.  Time-frequency  techniques  allow 
the  tracking  of  these  effects. 

1.2.2  Dynamic  Targets 

The  basic  Fourier  transform  assumes  signal  stationarity  in  the  sense  that  the  frequencies 
of  the  signal  components  do  not  change  with  time.  For  a  dynamic  target,  this  is  a  poor 
assumption.  Using  time-frequency  techniques  to  analyze  non-stationary  targets  maintains 
the  signal’s  relationship  to  time  while  displaying  how  the  frequency  content  of  the  signal 
varies  with  respect  to  time.  This  shows  how  the  target  movements  drive  the  scattered 
spectrum. 

1.2.3  Noise  Corrupted  Signals 

When  considering  signal  processing  techniques,  the  usefulness  of  these  techniques  has  to 
be  evaluated  with  respect  to  contaminated  signals.  Most  radar  and  remote  sensing  systems 
are  noise-limited  in  their  performance  [36].  It  was  important  to  examine  how  well  time- 
frequency  analysis  techniques  operate  in  a  noise-filled  environment.  These  distributions 
tend  to  be  very  sensitive  to  noise  when  trying  to  track  scattering  modes  containing  very 
low  signal  levels. 

1.2.4  Customized  Kernels 

In  most  sections  of  the  dissertation,  the  signals  were  processed  using  standard  kernels 
such  as  the  Wigner,  the  spectrogram,  and  the  binomial.  This  section  looks  at  the  use  of 
customized  kernels  to  generate  the  time-frequency  representation.  A  genetic  algorithm  was 
developed  which  maximized  the  third-order  Renyi  entropy  for  the  time-frequency  distribu¬ 
tion.  This  customized  kernel  provided  performance  similar  to  a  binomial  and  consistently 
exceeded  the  binomial  and  the  Wigner  distribution  when  evaluated. 


3 


1.3  Summary 


Time-frequency  techniques  clearly  demonstrated  their  value  for  analyzing  radar  based 
signals.  For  stationary  targets  examined  with  a  swept  frequency  signal,  these  techniques 
can  isolate  different  scattering  centers  and  scattering  modes.  For  dynamic  targets,  they 
provide  a  view  of  the  target’s  dynamics  unavailable  in  regular  Fourier  analysis. 
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CHAPTER  2 


BACKGROUND 


This  chapter  outlines  the  basic  mathematics  developed  and  used  in  this  research.  Start¬ 
ing  with  the  Fourier  transform  to  move  between  the  time  and  the  frequency  domains,  the 
idea  of  changing  domains  to  gain  insight  for  a  signal  is  extended  to  joint  time-frequency 
and  frequency-time  domains. 

The  chapter  also  examines  the  history  of  time-frequency  analysis  and  discusses  earlier 
developments  and  applications  of  these  techniques. 

2.1  Traditional  Spectral  Analysis 

The  Fourier  transform  and  inverse  Fourier  transform  are  a  set  of  relationships  which 
allow  representing  a  function  in  two  different  domains.  The  time,  or  temporal,  domain 
uses  time  as  the  independent  variable.  For  the  frequency,  or  spectral,  domain,  frequency 
serves  as  the  independent  variable.  Converting  a  time  function  into  its  frequency  domain 
representation  simplifies  many  operations  and  provides  insight  into  different  aspects  of  the 
signal  that  may  be  difiicult  to  discern  from  the  time  domain  representation. 

2.1.1  Standard  Usage 

The  Fourier  transform  is  given  by 

/OO 

5(T)e-^2-/-dr  (2.1) 

-OO 

The  corresponding  inverse  transform  is 

/OO 

G{X)e^^^^^dX  (2.2) 

-OO 

Different  disciplines  employ  various  permutations  of  the  Fourier  transform  relationships. 
By  using  /  to  represent  the  frequency  (in  hertz),  the  transform  equations  above  avoid  a 
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factor  of  l/27r  that  appears  when  working  with  uj  (in  radians/sec)  as  the  frequency  vari¬ 
able  [15]. 

For  discrete  signals,  the  relationship  takes  a  series  format  to  relate  the  frequency  domain 
series  G(p)  with  the  time  domain  series  g{n) 

0(9)  =  E  (2.3) 

n=0 

where  p  represents  the  frequency  index  and  n  represents  the  time  index. 

The  corresponding  inverse  transform  is 

»(')  =  vE0(9)e>"»  (2.4) 

^  g=0 

This  discrete  form  treats  the  original  g{n)  as  evenly  spaced  samples  of  a  periodic  signal. 
If  the  samples  are  taken  AT  seconds  apart,  and  there  are  N  samples  covering  one  period  of 
the  signal,  then  in  the  frequency  domain  the  samples  of  the  spectrum  will  be  spaced  1/AT 
hertz  apart,  giving  a  frequency  resolution  of  AT.  The  spectrum  will  also  be  periodic,  with 
N  distinct  values.  When  dealing  with  two-sided  spectra,  the  frequency  range  typically  is 
taken  from  -{N/2)AF  to  {N/2  -  1)AT. 

2.1.2  Power  Spectral  Density 

The  signal  information  in  the  frequency  domain  is  split  between  the  real  and  imaginary 
portions  of  the  spectrum.  It  is  often  more  convenient  to  work  with  the  power  spectral 
density,  a  real-valued  representation  equivalent  to  the  magnitude  squared  values  of  the 
Fourier  transform.  It  can  be  obtained  calculating  the  spectrum  first 

«.(/) = (2-5) 

or  by  taking  the  Fourier  transform  of  the  signal’s  autocorrelation  function 

5.(/)  =  T[R,{t)]  (2.6) 

where  the  autocorrelation  is  defined  as 

"  1-00 ""  0  ^  0  "  0 

2.1.3  Limitations 

The  Fourier  transform  has  proved  itself  to  be  a  powerful  and  useful  analysis  tool,  but  it 
does  have  limitations  when  applied  to  complicated  signals.  Fourier  transforms  are  poorly 
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suited  for  non-stationaxy  signals.  Although  no  information  is  lost  when  transforming  a 
signal  across  domains,  certain  aspects  and  characteristics  can  be  obscured  and  difficult 
to  interpret.  For  instance,  a  particularly  difficult  aspect  to  distinguish  with  traditional 
Fomier  transforms  are  changes  in  the  signal’s  frequency,  either  due  to  the  signal  starting 
and  stopping,  or  moving  from  one  frequency  to  another. 

2.2  Time- Frequency  Analysis 

Several  methods  exist  for  extending  the  basic  ideas  of  the  Fourier  transform  to  ade¬ 
quately  represent  non-stationary  signals.  Instead  of  completely  removing  the  dependence 
on  the  original  independent  variable,  it  is  carried  over  and  combined  with  a  second  inde¬ 
pendent  variable.  So  instead  of  the  transform 

A[g{x)]  =  G{y)  (2.8) 

the  new  transformation  gives 

B[g{x)]  =  G{x,y)  (2.9) 

There  have  been  many  different  transformations  developed  to  provide  this  joint  fre¬ 
quency  representation.  A  large  family  of  these,  collectively  known  as  Cohen’s  class,  share 
many  characteristics  and  cover  the  most  popular  time-frequency  transformations  [11,  12]. 
These  transforms  can  be  calculated  via  the  ambiguity  domain  or  via  the  autocorrelation 
domain. 

2.2.1  Ambiguity  Domain  Signal  Transformations 

The  time-frequency  transformations  used  in  this  study  are  members  of  Cohen’s  class. 
These  transformations  generate  a  two-dimensional  function  of  time  and  frequency  which 
represents  the  signal  energy  per  unit  time  per  unit  frequency  [11].  Members  of  this  class 
have  the  general  form  [7] 

^)  =  ^  T)g  +  I)  9*  (/^  -  |)  dfidrd^  (2.10) 

and  move  from  the  time- valued  function,  g{t),  to  a  function  representing  the  distribution 
of  the  original  signal  in  a  mixed  time-frequency  domain.  The  resulting  function  of  time 
and  frequency,  Cg{t,uj,(f>)  depends  on  <^(a;,  r),  the  transformation  kernel  chosen.  This  can 
be  modified  slightly  to  use  non-radian  frequencies  by  substituting  lj  — »  27r/  and  ^  2'kv. 
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Applying  this  to  the  previous  equation  gives 

A  Cohen’s  class  transformation  can  be  viewed  as  having  three  steps:  moving  from 
g{t)  — »  Ag{i/,T),  applying  a  weighting  fimction,  and  using  a  double  transformation  to  go 
from  Ag{i/,T)  —>■  Cg{t,f). 

The  first  step  moves  the  signal  from  the  time  domain  into  a  two-dimensional  domain, 
similar  to  the  ambiguity  domain  used  in  radar  signal  analysis  [36,  47].  In  the  radar  context, 
r  represents  a  time  shift,  and  /  is  a  frequency  shift. 

/OO 

g{t)g*{t  +  r)e^^^f^dt  (2.12) 

'OO 

The  first  sub- transfer  mat  ion  with  Cohen’s  class  gives  something  slightly  different  from 
the  definition  in  equation  2,12.  It  is  called  the  symmetric  ambiguity  function  and  is  defined 
to  be 

(*">  t)  =  5  +  0  9*  (9'  -  0  ^2.13) 

=  (2.14) 

Once  the  signal  is  moved  to  this  ambiguity  domain,  a  kernel  function,  ^(i/,  r),  is  ap¬ 
plied.  The  choice  of  this  kernel  determines  what  kind  of  time-frequency  transformation  is 
used.  Table  2.1  lists  several  of  the  distributions  discussed  by  Cohen  [12],  along  with  their 
transformation  kernels  and  some  properties. 


Distribution 

Name 

Kernel 

Resolution 

Time  Frequency 

Cross 

Terms 

Wigner-Ville 

1 

High 

High 

High 

Spectrogram 

J  h*{u  —  -f  §)dn 

Low 

High 

Low 

Exponential 

High 

High 

Low 

Table  2.1:  Selected  TF  Kernels 


Each  of  the  three  distributions  listed  has  desirable  properties  and  liabilities.  The  spec¬ 
trogram,  also  called  the  short  time  Fourier  transform,  involves  a  direct  tradeoff  between 
time  and  frequency  resolution.  Depending  on  the  chosen  window,  /i(t),  this  time-frequency 
distribution  can  provide  very  high  frequency  resolution  but  only  at  the  cost  of  poor  tempo¬ 
ral  resolution.  It  often  provides  an  unrealistic  viewpoint  of  the  time-frequency  structure  of 
complex  signals  [40]. 
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The  second  distribution,  the  Wigner,  was  introduced  in  the  context  of  quantum  me¬ 
chanics  [42]  and  later  adapted  and  applied  to  signal  processing  [39].  The  Wigner-Ville 
distribution  in  general  provides  the  highest  time  and  frequency  resolutions.  The  drawback 
when  using  this  form  is  the  presence  of  significant  cross  terms  between  all  possible  groups  of 
actual  signal  components.  This  causes  significant  clutter  and  confusion  when  complicated, 
multi-component  signals  are  analyzed  [13,  14,  19]. 

The  exponential  distribution  is  a  Reduced  Interference  Distributions,  or  RID.  The  goal  of 
a  RID  is  to  maintain  high  time  and  frequency  resolution  while  minimizing  the  contribution 
from  cross  terms  [7]. 

2.2.2  Autocorrelation  Domain  Signal  Transformations 

In  many  cases,  it  is  convenient  to  work  with  the  transformation  in  a  different  form. 
Instead  of  moving  the  signal  into  the  ambiguity  domain,  the  kernel  can  be  applied  using  a 
convolution  in  the  autocorrelation  domain.  In  this  case,  the  time  frequency  distribution  is 
given  by  [21] 


roo 

J—oo 

(2.15) 

using 

^  /  oo  ^  ^  ~  0  ~ 

(2.16) 

and 

^OO 

J  —  OO 

(2.17) 

By  defining  the  local  autocorrelation  function 

II 

+ 

* 

1 

(2.18) 

this  can  also  be  written  as 

poo 

Cg{t,  /;  V’)  =  /  Rg(t,  t) 

J—OO 

(2.19) 

Knowing  the  autocorrelation  domain  representation  of  the  kernel  allows  one  to  determine 
the  joint  time-frequency  representation  with  only  two  computationally  intensive  operations, 
a  convolution  and  a  Fourier  transform,  instead  of  the  three  Fourier  transforms  required  by 
the  original  formulation  in  equation  2.11. 

2.2.3  Sampled  Signal  Analysis 

Since  most  work  involves  sampled  data,  the  transformations  take  on  slightly  different 
discrete  forms,  with  discrete  transforms,  as  shown  in  equations  2.3  and  2.4,  replacing  the 


9 


continuous  transforms  in  equations  2.11  and  2.19.  The  ambiguity  domain  based  expression 
becomes 


.  P  N  N  , 

=  Jr  EEEs  ("i+  -)  g*  ("i-  2)  ^(^’0 

i=0g=0m=0  \  \ 


^3^qm 


(2.20) 


while  the  autocorrelation  domain  expression  is  given  by 


p-i 


Cg{n,p) 


i=0  fc--oo 


(2.21) 


even  I 


These  formulations  exhibit  the  same  basic  properties  as  seen  in  their  continuous  coun¬ 
terparts  [21].  A  major  difference  is  the  sampled  signal  is  assumed  periodic,  with  the  period 
length  matching  the  signal  vector  g{n).  The  overall  effect  is  to  generate  a  two-dimensional 
time-frequency  distribution  which  is  periodic  in  both  directions.  Properly  defining  the 
summation  limits  is  essential  to  avoid  problems  of  aliasing. 


2.3  Desirable  Time-Frequency  Distribution  Properties 

When  considering  different  time-frequency  distributions,  there  are  several  desirable 
properties  we  would  like  the  distribution  to  exhibit.  They  range  in  importance,  with  dif¬ 
ferent  kernels  satisfying  different  properties.  Ideally,  this  joint  distribution  would  behave 
as  a  power  or  energy  density  function  and  would  share  properties  found  in  two-dimensional 
probability  density  functions.  Table  2.2  lists  these  properties. 

2.3.1  Non-Negative 

Since  the  time-frequency  distribution  can  be  viewed  as  an  energy  density  distribution, 
it  does  not  make  sense  for  any  values  in  the  time-frequency  plane  to  be  negative.  If  these 
negative  values  can  be  accepted  as  a  numerical  artifact,  loosening  this  requirement  makes 
meeting  the  other  properties  possible. 

The  only  common  time-frequency  kernel  to  support  this  property  is  the  spectrogram  [22]. 

2.3.2  Real  Valued 

Since  the  time-frequency  distribution  is  related  to  the  energy  spectral  density,  all  values 
should  be  real.  To  achieve  this,  the  time-frequency  kernel  possesses  Hermitian  symmetry 
with  respect  to  v  and  r  in  the  ambiguity  domain 

<j>{u,T)  =  <t>{-u,-T)  (2.22) 
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PO  Non-Negative 

PI  Real  Valued 

P2  Time  Shift  Invariant 

P3  Frequency  Shift  Invariant 
P4  Time  Marginal 

P5  Frequency  Marginal 

P6  Instantaneous  Frequency 
P7  Group  Delay 

P8  Time  Support 

P9  Frequency  Support 

PIO  Reduced  Interference 

Table  2.2:  Desirable  Distribution  Properties 

and  possesses  Hermitian  symmetry  with  respect  to  r  in  the  autocorrelation  domain 

(2.23) 


2.3.3  Time  Shift  Invariant 

If  the  original  signal  is  shifted  in  time,  the  time-frequency  distribution  should  exhibit 
an  equivalent  shift  along  the  time  axis.  This  condition  will  be  met  with  any  kernel  in  a 
Cohen’s  class  transformation. 

2.3.4  Frequency  Shift  Invariant 

If  the  original  signal  is  shifted  in  frequency,  the  time- frequency  distribution  should  ex¬ 
hibit  an  equivalent  shift  along  the  frequency  axis.  The  Cohen’s  class  transformations  auto¬ 
matically  satisfy  this  condition  regardless  of  the  kernel  selected. 

2.3.5  Time  Marginal 

Paralleling  the  concepts  of  a  two-dimensional  probability  density  function,  the  time- 
frequency  distribution  should  give  the  correct  value  for  the  instantaneous  signal  power  at 
any  time  by  integrating  across  the  frequency  axis 

/oo 

Cg{tJ)df  (2.24) 

-OO 

For  this  condition  to  be  met,  the  ambiguity  kernel  must  be  equal  to  1  along  the  frequency 
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axis,  and  the  autocorrelation  kernel  must  be  a  delta  function  along  the  time  axis. 


=  1 


=  6{t) 


(2.25) 

(2.26) 


2.3.6  Frequency  Marginal 


As  just  mentioned  for  the  time  marginal,  the  power  spectral  density  should  result  after 
integrating  out  the  dependence  on  t. 


/OO 

C,(t,f)dt 

-OO 


(2.27) 


This  condition  can  be  met  using  ambiguity  kernels  having  the  value  of  1  along  the  lag 
axis  {u  =  0)  or  equivalently,  an  autocorrelation  domain  kernel  which  when  integrated  with 
respect  to  t  will  equal  1  for  all  values  of  t. 


(p{0,  t) 


/OO 

V’(t,  T)dT  = 

-OO 


(2.28) 

(2.29) 


2.3.7  Instantaneous  Frequency 


The  average  frequency  at  any  time  should  equal  the  instantaneous  frequency.  Treating 
the  time-frequency  distribution  like  a  probability  distribution  allows  expressing  the  average 
frequency  as 


favg{^)  — 


S-^CM,f)df 


and  the  instantaneous  frequency 


(2.30) 


(2.31) 


where  0{t)  is  the  signal  phase  at  time  t. 

To  achieve  this  property,  the  kernel  must  satisfy  several  conditions.  It  must  provide  a 
proper  time  marginal  (equations  2.25  or  2.26).  Also,  the  partial  derivative  of  the  ambiguity 
domain  kernel  taken  with  respect  to  the  lag  variable,  r,  and  evaluated  at  r  =  0  must  be 


zero  for  all  u 


d<t>{v,T) 

dr 


(2.32) 
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2.3.8  Group  Delay 

For  a  given  frequency,  the  average  time  across  the  time-frequency  distribution  should 
equal  the  group  delay  of  the  original  signal.  The  average  time  is  given  by 


tavgif)  — 


and  the  group  delay 


tgif)  = 


1  de{f) 


(2.33) 


(2.34) 


27r  df 

with  9{f)  giving  the  signal  phase  as  a  function  of  frequency. 

To  achieve  this  property,  the  kernel  must  satisfy  several  conditions.  It  must  provide 
a  proper  frequency  marginal  (equations  2.28  or  2.29).  Also,  the  partial  derivative  of  the 
ambiguity  domain  kernel  taken  with  respect  to  the  frequency  shift  variable,  u,  and  evaluated 
at  =  0  must  be  zero  for  all  r 

dHi',  t)  I 


du 


=  0 


j/=0 


(2.35) 


2.3.9  Time  Support 


The  time  support  property  states  if  a  signal  has  non-zero  values  only  over  the  range 
|t|  <  tc,  then  the  time-frequency  distribution  should  only  have  non-zero  values  over  the 
same  range. 

This  property  corresponds  to  an  ambiguity  domain  kernel  satisfying 


/OO 

-OO 


or  an  autocorrelation  kernel  where 


for  |t|  <  2|t|. 


ip{t,  r)  =  0 


(2.36) 


(2.37) 


2.3.10  Frequency  Support 

Frequency  support  means  if  the  spectrum  G{f)  is  zero  for  all  |/|  >  /c,  then  the  time- 
frequency  distribution  Cg{t,f)  must  also  be  zero  when  |/|  >  fc- 

The  kernel  required  for  this  property  has  an  ambiguity  domain  kernel  where 


/OO 

-OO 


(2.38) 


for  all  \v\  <  2|/|. 
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2.3*11  Reduced  Interference 


The  final  desirable  property  for  a  time-frequency  distribution  is  to  have  small  inter¬ 
ference  terms  between  the  signal  components.  While  not  precisely  defined  as  the  previous 
properties,  this  is  very  important  since  interference  terms  can  seriously  hinder  interpretation 
of  the  time-frequency  distribution. 

To  get  reduced  interference  terms,  the  ambiguity  domain  representation  of  the  kernel 
must  behave  as  a  low-pass  filter  with  respect  to  both  the  frequency  shift,  and  the  time 
lag,  r,  axes. 


2.4  Kernel  Design  Procedures 


Although  the  properties  listed  in  table  2.2  look  formidable,  it  is  possible  to  meet  all  the 
properties  except  PO  (non-negativity)  through  a  straightforward  design  procedure  [22]. 

The  basic  steps  involve  taking  a  smooth  elementary  function,  h{t),  defined  over  the 
region  —l/2<t<  1/2.  The  function  should  enclose  an  area  of  1,  be  symmetric  about 
t  =  0,  and  go  smoothly  to  zero  as  \t\  approaches  1/2. 

Once  this  elementary  function  is  defined,  the  ambiguity  domain  kernel  is  given  by  the 
Fourier  transform,  with  the  product  ur  substituted  for  the  frequency  variable 


0(1/,  r)=  J 


f~i/r 


(2.39) 


Moving  back  into  the  autocorrelation  domain,  using  the  properties  of  Fourier  and  inverse 
Fourier  transforms,  the  kernel  becomes 


1 


(2.40) 

(2.41) 


2.5  Alias-Free  Formulation 

When  evaluating  a  discrete  autocorrelation  function,  a  problem  arises  due  to  the  ar¬ 
guments  m  +  1/2  and  m  —  1/2.  While  they  cause  no  problems  for  continuous  signals,  a 
discrete  signal  is  typically  treated  as  undefined  (or  zero)  for  non-integer  arguments.  This 
limitation  has  the  effect  of  zeroing  out  the  autocorrelation  for  all  odd  lag  values.  When 
viewed  along  the  lag  axis  (/),  instead  of  having  P  points  separated  by  AT,  there  are  only 
P/2  points  separated  by  2AT.  The  end  result  is  undersampling  the  distribution  by  a  factor 
of  2,  allowing  a  frequency  span  of  only  —{N/A)AF  to  (iV/4  —  1)AF  before  aliasing  occurs. 
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To  avoid  this  problem,  it  is  possible  to  re-formulate  the  definition  of  the  autocorrelation 
function  to  provide  a  half  step  shift  to  place  the  values  back  on  the  sample  points.  Phasing 
problems  are  avoided  by  incorporating  a  corresponding  negative  half  step  into  the  kernel 
definition. 

2.6  Dual  of  Transformations 

A  quick  examination  of  equations  2.1  and  2.2  shows  the  forward  and  inverse  Fourier 
transforms  differ  only  in  the  sign  of  the  exponential  term.  This  means  any  properties  or 
relationships  exhibited  by  the  forward  transform  will  also  apply  to  the  inverse  transform 
simply  by  substituting  g{—t)  for  G{f). 

With  the  discrete  transforms  given  by  equations  2.3  and  2.4,  the  same  relationship  almost 
holds.  The  only  difference  is  a  scaling  factor  1/N  that  applies  to  the  inverse  transform. 

This  similarity  between  forward  and  inverse  transforms,  or  duality  [15],  allows  recasting 
the  time-frequency  transformations  into  frequency-time  transformations.  Given  a  function 
of  frequency,  one  may  find  the  temporal  response  as  a  function  of  frequency.  Problems  of 
this  type  are  typical  in  measured  electromagnetic  scattering  data. 

2.7  Emphasized  Kernels 

As  presented,  Cohen’s  class  covers  a  broad  range  of  time-frequency  distributions.  Em¬ 
phasis  has  been  placed  on  four  different  distributions:  the  spectrogram,  the  Wigner  distri¬ 
bution,  the  binomial  distribution,  and  a  genetically  derived  distribution. 

2.7.1  Spectrogram 

The  spectrogram,  also  known  as  the  short  time  Fourier  transform,  or  STFT,  provides  the 
most  intuitive  time  frequency  distribution.  While  it  can  be  cast  in  the  standard  Cohen’s 
class  form  with  the  kernel  shown  in  table  2.1,  the  usual  approach  centers  a  window  on 
the  time  of  interest  and  takes  the  Fourier  transform  of  the  windowed  data.  Finding  the 
magnitude  squared  of  the  result  completes  the  transformation. 

The  resolution  provided  by  this  transform  depends  on  the  length  of  the  window.  A  long 
window  provides  good  frequency  resolution,  but  it  degrades  the  transform’s  ability  to  track 
the  temporal  position  of  the  signal.  A  narrow  window  improves  the  ability  to  determine  the 
time  of  a  signal  component  accurately,  but  this  comes  at  the  expense  of  coarse  frequency 
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resolution. 


The  spectrogram  is  best  suited  for  signals  with  a  slowly  varying  frequency  content. 
While  the  spectrogram  provides  a  non-negative  time-frequency  representation,  it  does  not 
satisfy  the  time  or  frequency  marginals. 


2.7.2  Wigner  Distribution 


The  Wigner  distribution,  also  referred  to  as  the  Wigner-Ville  distribution,  was  originally 
used  in  the  study  of  quantum  statistical  mechanics,  working  with  position  and  momentum 
instead  of  time  and  frequency  [42].  Later  work  by  Ville  independently  derived  this  trans¬ 
formation  for  signal  processing,  using  time  and  frequency  [39].  The  ambiguity  domain 
representation  for  the  Wigner  is  a  delta  function  along  the  time  axis,  t,  and  is  independent 
of  the  time  lag,  r. 

The  delta  functional  form  of  the  Wigner  generates  two  properties  for  the  distribution. 
Having  a  delta  function  in  the  kernel  makes  it  unnecessary  to  carry  out  the  convolution 
with  the  kernel,  making  the  transform  quick  to  calculate.  This  also  provides  the  highest 
possible  resolution  in  the  time  frequency  domain. 

The  Wigner  distribution  possesses  nine  of  the  eleven  desirable  properties  listed  in  ta¬ 
ble  2.2.  The  only  two  it  fails  to  meet  are  non-negativity,  PO,  and  it  is  not  a  reduced 
interference  distribution,  PIO.  A  drawback  of  the  Wigner  distribution  is  the  lack  of  any 
smoothing  in  the  time-frequency  domain.  Large  cross  terms  exist  even  when  the  true  com¬ 
ponents  are  clearly  separated.  While  this  is  tolerable  when  few  signal  components  exist, 
the  number  of  cross  terms,  Nx,  generated  from  N  actual  components  is  given  by 


_  N{N-1) 


(2.42) 


2.7.3  Binomial  Distribution 

The  binomial  distribution  is  typically  generated  using  a  discrete  autocorrelation  domain 
kernel  which  is  closely  related  to  the  continuous  exponential  distribution  in  table  2.1.  The 
binomial  distribution  sacrifices  some  of  the  resolution  provided  by  the  Wigner  to  reduce 
the  cross  terms,  similar  to  how  a  tapered  anterma  array  trades  off  pattern  beamwidth  for 
reduced  sidelobes. 

The  binomial  kernel  gives  a  reduced  interference  distribution  which  satisfies  all  the 
desired  properties  except  non-negativity. 
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2.7.4  Genetic  Kernel 


Chapter  7  discusses  a  genetic  algorithm  used  to  design  a  time-frequency  distribution. 
The  algorithm  selected  coefficients  for  a  sinusoidal  representation  of  the  elementary  function 
h{t)  as  outlined  in  section  2.4. 

The  distributions  found  via  this  method  were  very  similar  to  those  provided  by  the 
binomial.  Like  the  binomial,  they  met  all  the  desirable  properties  except  non-negativity. 
This  kernel  exceeded  the  binomial  performance  when  evaluating  the  entropy  of  the  resulting 
distribution.  Chapter  7  gives  more  information  on  this  kernel. 

2.8  Review  of  Time-Prequency  Developments 

The  basics  of  time-frequency  techniques  trace  back  to  Fourier  analysis  through  the  short 
time  Fourier  transform.  The  spectrogram  began  seeing  extensive  use  for  speech  signals  in 
the  1940’s  [8,  31]  due  to  the  inadequacy  of  simple  Fourier  transforms  for  analyzing  non¬ 
stationary  signals. 

In  1932,  Wigner  introduced  a  means  of  estimating  momentum  and  position  in  quantum 
mechanics  [42].  The  basic  techniques  were  applied  in  a  signal  processing  context  by  Ville  in 
1948  [39],  leading  this  approach  to  often  be  referred  to  as  the  Wigner-Ville  distribution  [32]. 
This  had  the  advantage  of  much  higher  resolution  than  the  spectrogram. 

Alternative  representations  for  the  joint  time-frequency  domain  were  developed  through 
the  1950’s  and  1960’s  with  researchers  such  as  Page,  Levin,  Rhihaczek,  and  Born-Jordan. 
The  widely  varying  characteristics  of  these  different  distributions  were  consolidated  in  1966 
by  Cohen  into  a  general  class  of  transformation  for  joint  distributions  [11].  With  this 
representation,  an  arbitrary  kernel  (j){0,T)  differentiated  the  different  distributions.  This 
formulation  allowed  analysis  and  comparisons  between  kernels  and  led  to  experimentation 
on  new  kernels  [12]. 

In  1980,  Claasen  and  Mecklenbrauker  brought  the  Wigner  distribution  to  the  forefront 
with  a  series  of  articles  discussing  the  strengths,  weaknesses,  and  applications  for  this  anal¬ 
ysis  approach  [8,  9,  10]. 

A  lot  of  recent  work  has  led  to  a  better  understanding  of  the  fundamental  limits  of 
performance  for  Cohen’s  class  distribution  and  finding  kernel  and  computational  techniques 
to  provide  the  clearest  time-frequency  representations  [2,  18,  22,  21,  45]. 
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2.9  Time-Frequency  Applications 


Time-frequency  representations  have  been  used  extensively  for  analyzing  non-stationary 
systems  where  the  signal  spectra  varies  as  a  function  of  time.  In  some  cases,  the  objective 
was  to  characterize  individual  components  of  the  signal.  In  other  cases,  the  overall  patterns 
in  the  time-frequency  representation  served  as  the  metric  for  making  decisions  about  the 
signal. 

2.9.1  Marine  Mammal 

A  large  amoimt  of  work  has  been  conducted  analyzing  the  whistles  and  other  sounds 
generated  by  marine  mammals.  Decomposing  these  signals  in  the  time-frequency  domain 
shows  dolphins  and  whales  using  frequency  modulated  sonar  signals  employing  sophisticated 
techniques  [1]. 

2.9.2  Biomedical 

Another  area  where  time-frequency  techniques  have  been  applied  heavily  is  the  analysis 
of  other  biological  signals.  Some  of  these  signals  have  been  acoustic,  such  as  differenti¬ 
ating  benign  jaw  clicks  from  those  requiring  corrective  treatment  [1].  Other  medical  uses 
have  examined  the  body’s  own  electrical  signal,  such  as  trying  to  characterize  the  electroen¬ 
cephalogram  (EEG)  signals  of  brainwave  activities  associated  with  epileptic  seizures  [1].  The 
heart’s  electrical  signals,  measured  with  an  electrocardiogram  (EKG)  have  been  analyzed 
to  try  and  differentiate  the  signals  of  a  healthy  heart  from  a  diseased  heart  [43]. 

2.9.3  Mechanical  Engineering 

Time-frequency  techniques  are  also  finding  uses  in  disciplines  such  as  mechanical  en¬ 
gineering.  Researchers  have  monitored  vibrations  from  cutting  tools,  trying  to  determine 
characteristics  in  the  tool’s  chatter  that  may  indicate  imminent  failure  [1]. 

2.9.4  Electromagnetics 

Cohen’s  class  of  time-frequency  distributions  have  not  been  used  heavily  in  the  field  of 
electromagnetics.  Some  work  has  employed  spectrograms  and  wavelet  analysis  [23,  25,  28], 
but  high  resolution  time-frequency  techniques  like  the  Wigner  and  the  binomial  have  been 
largely  ignored.  The  typical  application  has  been  the  analysis  of  scattering  modes  when 
analyzing  a  target  illuminated  by  a  wideband  system. 
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2.10  Summary 


Time-frequency  analysis  has  been  largely  ignored  in  the  study  of  electromagnetic  scat¬ 
tering.  While  some  related  work  has  been  done,  the  emphasis  has  been  on  using  the  wavelet 
transform  instead  of  Cohen’s  class  distributions  such  as  the  Wigner  or  binomial  distribu¬ 
tions.  The  typically  frequency  dependent  nature  of  scattering  and  the  time-varying  spec¬ 
tra  associated  with  dynamic  targets  make  time-frequency  techniques  an  attractive  tool  for 
working  with  scattering  problems. 
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CHAPTER  3 


STATIONARY  TARGETS 


3.1  Introduction 

Anechoic  chambers  can  be  used  to  measure  both  the  spatial  pattern  and  the  frequency 
response  of  a  target.  As  long  as  the  measurement  system  provides  a  stable,  adjustable 
frequency  source,  the  chamber  can  also  provide  the  information  necessary  to  evaluate  the 
impulse  response  of  the  target  and  chamber. 

3.2  Background 

Anechoic  chamber  measurements  to  determine  the  pulse  response  of  a  target  normally 
do  not  use  an  actual  pulse.  Rather  than  deal  with  the  physical  difficulties  of  generating  high 
power,  short  duration  pulses,  a  measurement  system  illuminates  the  target  with  a  stepped 
continuous  wave  (CW)  signal.  Instead  of  trying  to  accurately  record  and  analyze  a  pulse 
return  in  real  time,  the  frequency  is  selected,  the  system  is  allowed  to  reach  steady  state, 
and  a  complex-valued  data  point  (magnitude  and  phase),  G{fn)  is  collected  before  moving 
on  to  the  next  frequency  [24,  41]. 

Using  the  inverse  Fourier  transform  gives  the  target  and  chamber  response  as  a  function 
of  time 

9itn)=T-^[G{fn)]  (3.1) 

which  approximates  the  impulse  response. 

If  N  samples  are  taken,  each  separated  by  A/  hertz,  these  samples  can  be  transformed 
to  the  time  domain  where  they  will  provide  N  samples  of  the  impulse  response  over  a  total 
temporal  span  of  1/A/  seconds. 

Assuming  the  wave  propagates  at  a  velocity  equal  to  the  speed  of  light  in  free  space. 
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t;  =  c  =  3  X  10®Tn/s.  This  means  the  time  response  can  be  expressed  as  a  range  response 


9{r)  =  (3.2) 

This  same  approach  to  finding  the  time  response  can  also  be  applied  to  numerically 
simulated  results.  Finding  the  spectral  response  provides  the  sampled  impulse  response  of 
a  virtual  anechoic  chamber  with  a  maximum  imambiguous  range 


RmitK  — 


Measurements  within  this  virtual  chamber  have  a  range  resolution  determined  by  the  maxi¬ 
mum  and  minimum  frequencies  used  in  calculating  the  transform.  For  N  frequencies  spaced 
A/  apart,  the  range  resolution  is 


Ar  = 


V 

2(iV-l)A/ 


(3.4) 


3.3  Duality  and  Time-Frequency  Transformations 

The  typical  time-frequency  transformations  move  a  function  of  time  into  a  time-frequen¬ 
cy  space  using  the  Cohen’s  class  transformation 

Cg{t,  f,  ci>)  =  /_"  r)g  -h  0  g*  (u  -  0  dudrdv  (3.5) 

This  can  also  be  viewed  as  a  Fourier  transform  to  convert  the  signal’s  local  autocorrela¬ 
tion  into  an  ambiguity  function,  applying  a  weighting  function,  or  kernel,  then  doing  a  two- 
dimensional  Fourier  transform  to  move  from  the  ambiguity  domain  to  the  time-frequency 
domain. 

Cg{tJ)  =  Trf  [rzlt  [5  («  +  0  9*  (^  -  0]  <^('^-^)})  (3-6) 

Typically,  instead  of  transforming  the  local  autocorrelation  into  the  ambiguity  domain, 
the  kernel  is  applied  by  performing  a  convolution  along  the  time  axis,  t,  in  the  autocorre¬ 
lation  domain.  After  convolving,  taking  the  Fomier  transform  of  the  lag  variable,  r,  yields 
the  time  frequency  distribution. 

Cg{t, /)  =  Trf  [Rg{t, t)  (8>t  V’(t,  r)]  (3.7) 

with  Rg{t,T)  as  the  local  autocorrelation  and  <g)t  representing  the  convolution  with  respect 
to  t. 

When  beginning  with  a  frequency  domain  signal,  G{f),  the  spectral  autocorrelation  is 

flG(t',/)  =  G(/  +  0G'(/-5)  (3.8) 
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Taking  the  inverse  Fourier  transform  with  respect  to  /  converts  the  spectral  autocorre¬ 
lation  into  the  same  ambiguity  domain  representation  that  would  have  resulted  if  starting 
with  g{t).  This  means  the  time-frequency  representation  for  G{f)  is  given  by 

Ca{t,f)  =  [g  (#.+  0  G‘  (p  -  0]  (3.9) 


This  result  can  also  be  obtained  by  applying  an  inverse  Fourier  transform  on  the  fre¬ 
quency  shift  variable,  i/,  from  the  spectral  autocorrelation  after  convolving  with  the  spectral 
autocorrelation  kernel,  (r',/) 


Cait,  f)  =  [Rair^,  f)  /)]  (3.10) 


Working  through  these  transformations,  the  duality  between  forward  and  inverse  Fourier 
transforms  provides  a  simplification.  Substituting  G{f)  for  g{t)  and  working  through  the 
temporal  time-frequency  distribution  steps  will  generate  the  result  Cg  (/,<).  This  simple 
switching  of  axis  works  provided 

=  (3.11) 


This  rotational  symmetry  will  happen  whenever  the  ambiguity  domain  kernel  is  a  product 
kernel  such  that 

(j>{u,  t)  -  (t>{vT)  (3.12) 

All  the  kernels  used  in  this  chapter  satisfy  that  condition. 


3.4  Scattering  and  Dispersion 

When  considering  the  scattered  radar  retmrn  from  a  target,  a  variety  of  modes  may 
contribute.  Based  on  the  target  composition  and  physical  arrangement  it  is  possible  to 
separate  the  returns  from  scattering  centers  which  are  in  different  range  bins  on  the  target. 

Scattering  modes  may  be  divided  into  dispersive  and  non-dispersive  modes.  A  non- 
dispersive  mode  will  behave  identically  at  all  frequencies.  A  dispersive  mode  will  depend 
on  a  mechanism  which  is  frequency  dependent  [33]. 

When  processed  using  a  standard  Fourier  transform,  sharp  resolution  is  possible  on  non- 
dispersive  modes.  The  energy  returned  through  the  dispersive  modes  will  appear  spread 
across  a  band  of  range  bins  since  the  signal  travels  at  a  non-constant  velocity  which  changes 
as  a  function  of  frequency.  However,  using  a  frequency-time  transformation  clearly  shows 
the  difference  between  the  dispersive  and  non-dispersive  modes,  and  can  also  differentiate 
various  dispersive  modes. 
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3.4.1  Reflective  Returns 


Reflective  returns  follow  from  geometric  optics  and  are  found  by  tracing  simple  reflec¬ 
tions  off  the  target  surface  with  the  angle  of  reflection  being  equivalent  to  the  angle  of 
incidence  relative  to  the  local  normal  on  the  target.  These  returns  can  involve  single  reflec¬ 
tions  directly  back  to  the  radar  or  multiple  bounce  reflections.  Over  the  range  of  frequencies 
considered  these  modes  are  non-dispersive. 

3.4.2  Creeping  Waves 

Another  scattering  mode  considered  involves  creeping  waves  which  attach  themselves  to 
the  target  and  travel  around  the  target  before  shedding  energy  back  toward  the  radar.  As 
these  waves  travel  around  the  target,  the  velocity  of  propagation  as  well  as  the  energy  shed 
toward  the  radar  will  vary  as  a  function  of  frequency.  At  low  frequencies,  the  waves  travel 
closely  attached  to  the  surface,  while  at  higher  frequencies  the  waves  are  offset  slightly  from 
the  surface.  Using  a  frequency-time  transform,  the  position  of  these  scattering  centers  will 
move  with  frequency.  The  slope  will  be  determined  by  how  far  the  energy  travels  while  in 
a  dispersive  mode. 

These  modes  are  difficult  to  track  accurately  because  they  carry  much  less  energy  than 
the  specular  returns  for  simple  targets. 

3.4.3  Combined  Modes 

Finally,  there  are  combined  modes  to  consider.  These  modes  have  energy  travel  a  path 
that  is  sometimes  specular  and  sometimes  dispersive.  They  can  also  be  categorized  by  how 
much  of  the  energy  path  is  along  a  dispersive  mode. 

3.5  Transforms  Considered 

For  this  study,  three  different  kernels  were  used  to  perform  the  frequency-time  transfor¬ 
mation.  The  kernels  were  chosen  because  of  their  widespread  use,  ease  of  implementation, 
and  the  variety  in  their  distribution  characteristics. 

3.5.1  Spectrogram 

This  transform  is  also  referred  to  as  the  short  time  Fourier  transform,  or  STFT.  It 
involves  moving  a  window  in  time  across  the  signal  and  taking  the  transform  for  each 
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window.  This  transform  gives  a  non-negative  distribution  with  good  noise  immunity  and 
low  cross  terms.  The  distribution  satisfies  neither  marginal  characteristic  and  generally 
provides  poor  resolution. 

3.5.2  Wigner 

While  the  spectrogram  is  one  of  the  most  intuitive  time-frequency  distributions,  the 
Wigner  employs  the  simplest  kernel,  a  delta  fimction.  The  result  is  a  very  high  resolution 
technique  which  suffers  from  high  interference  terms. 

3.5.3  Binomial 

The  binomial  distribution  was  developed  as  a  general  purpose  kernel  to  provide  high 
resolution  with  low  interference  terms.  While  the  binomial  cannot  match  the  resolution 
of  the  Wigner,  the  lack  of  large  interference  cross  terms  typically  results  in  a  clearer  joint 
domain  representation,  especially  with  complicated  signals 

3.6  Target  Geometries 

For  this  study,  the  basic  target  was  a  15  centimeter  radius  perfectly  conducting  sphere. 
It  was  chosen  due  to  the  canonical  nature  of  the  target.  It  was  simple  to  numerically  model, 
and  the  scattering  from  the  sphere  is  relatively  simple  and  well  understood.  The  sphere 
scatters  energy  using  two  basic  mechanisms:  a  simple  reflection  governed  by  geometrical 
optics,  and  dispersive  scattering  based  on  creeping  waves  traveling  around  the  surface  of 
the  sphere. 

3.6.1  Single  Sphere 

The  simplest  case,  shown  in  figure  3.1,  consists  of  a  single  perfectly  conducting  sphere 
in  free  space  illuminated  by  an  incident  plane  wave.  The  sphere  modeled  had  a  radius,  r, 
of  15  centimeters.  The  solution  for  the  scattered  field  was  obtained  by  evaluating  the  Mie 
series  [46]  at  the  frequencies  of  interest.  The  expected  returns  consist  of  a  reflective  return 
off  the  front  of  the  sphere  which  appears  a  distance  r  ahead  of  the  phase  reference  at  the 
sphere’s  center.  Most  of  the  reflected  energy  comes  from  this  mode.  The  next  contributor  is 
the  first  creeping  wave  which  travels  an  additional  irr  beyond  the  center  of  the  sphere.  Since 
the  range  values  are  based  on  two-way  propagation  times,  this  scattering  mode  appears  7rr/2 
behind  the  sphere  center.  The  other  detectable  mode  involves  a  second  creeping  wave,  one 
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Figure  3.1:  Single  Sphere 


which  travels  one-and-a-half  times  around  the  sphere  before  shedding  energy  back  towards 
the  radar.  Since  the  traveling  wave  continually  sheds  energy  as  it  travels  around  the  sphere, 
the  actual  impact  of  this  component  is  negligible. 

3.6.2  Two  Spheres:  In-line 

The  second  geometry,  shown  in  figure  3.2  involved  two  spheres,  each  with  a  radius,  r. 


1  meter 

4 - ► 


15  cm 

Figure  3.2:  Two  Spheres:  In-Line 


of  15  centimeters,  placed  with  a  distance,  d,  of  1  meter  between  their  centers.  Adding  the 
interactions  between  the  spheres  leads  to  several  more  modes  making  contributions  to  the 
scattered  signal.  The  primary  scattering  centers  are  the  front  face  of  each  sphere,  located 
r  in  front  of  the  phase  center  (the  center  of  the  leading  sphere)  and  d  —  r  behind  the  phase 
center.  The  next  modes  involve  the  first  creeping  wave  for  each  sphere.  These  appear  7rr/2 
behind  each  sphere  center,  or  at  ranges  of  7rr/2  and  d  —  7rr/2  behind  the  phase  center. 

When  examining  the  returned  signal,  the  scattering  modes  which  involve  creeping  waves 
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will  not  have  a  stationary  range  in  the  frequency-time  distribution.  The  energy  in  these 
dispersive  modes  will  also  drop  off  as  the  frequency  increases. 

3.6.3  Two  Spheres:  Broadside 

In  the  third  sphere  setup,  the  line  between  the  spheres  is  oriented  perpendicular  to 
the  direction  of  wave  propagation.  For  this  geometry,  the  polarization  of  the  illuminating 
wave  becomes  significant  when  considering  the  interactions  between  the  spheres  due  to  the 
behavior  of  creeping  waves. 

Figure  3.3  shows  the  two  spheres  at  a  broadside  orientation  to  the  source.  As  seen  in 
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Figure  3.3:  Two  Spheres:  Broadside 


the  previous  two  models,  the  specular  reflection  off  the  front  edge  of  the  sphere  will  be  the 
first  scattering  mode.  In  this  case,  the  energy  off  both  spheres  will  reinforce  each  other  and 
provide  a  single  large  return  a  distance  r  in  front  of  the  phase  reference.  The  combining 
effect  will  also  give  a  large  creeping  wave  return.  Both  these  modes  are  essentially  single 
sphere  modes  happening  independently  with  each  sphere.  They  do  not  depend  on  the 
incident  polarization. 

The  closest  scattering  mode  to  involve  interactions  between  the  spheres  is  the  reflection 
of  energy  off  one  sphere  towards  the  second,  then  off  the  second  back  in  the  direction  of  the 
somce.  This  mode  has  an  apparent  position  (d/2)  —  r  behind  the  phase  reference. 

Due  to  the  sensitivity  of  some  modes  to  the  polarization  of  the  illuminating  field,  two 
trials  were  computed.  The  first  trial  oriented  the  E-field  parallel  to  the  line  connecting 
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the  spheres.  This  allowed  sphere  interactions  based  on  creeping  waves.  The  second  trial 
oriented  the  El-field  perpendicular  to  the  line  between  the  spheres.  This  essentially  drove 
the  creeping  waves’  contributions  to  cross-sphere  effects  to  zero. 


3.7  Simulation  Parameters 


The  targets  outlined  above  were  simulated  numerically  using  two  different  packages:  a 
Mie  series  code  for  the  single  sphere  and  a  body  of  rotation  code  for  the  two-sphere  cases. 
In  each  case,  the  monostatic  backscatter  measured  included  the  magnitude  and  phase  of 
the  electric  field. 

For  the  two-sphere  targets,  the  spectrum  simulated  covered  the  band  from  1.00  gigahertz 
through  7.55  gigahertz  in  50  megahertz  steps.  This  provided  132  samples  for  processing. 
The  50  megahertz  sampling  interval  provided  a  maximum  unambiguous  range  scale  of  3 
meters  using  the  relationships 
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(3.13) 
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50  MHz 
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=  20  nsec 
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—  •  T 

2  -‘■max 


=  1.5  X  10®  —  •  20  nsec 
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(3.15) 


(3.16) 
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=  3  m 


(3.18) 


The  range  resolution,  Ar,  depends  on  the  span  of  frequencies  covered  by  the  spectral 
data.  With  spectral  information  covering  the  range  /min  to  /max?  the  inverse  Fourier  trans¬ 
form  can  provide  range  resolution  based  on 


Ar  = 


c 


2  (/max  /min) 


(3.19) 


assuming  the  signal  propagates  at  the  speed  of  light. 
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For  the  two-sphere  targets,  when  /max  ==  7.55  GHz  and  /min  =  1-00  GHz,  the  spatial 
resolution,  Ar,  is  2.29  centimeters. 

In  general,  the  resolution  in  the  frequency-time  domain  will  depend  on  the  length  of  the 
transform  used  after  applying  the  selected  kernel  in  the  autocorrelation  domain.  With  a 
frequency  sampling  interval  of  A/,  an  N  point  transform  equates  to 


Ar  = 


2{N  -  1)A/ 


(3.20) 


The  single  sphere  was  simulated  over  a  band  of  frequencies  starting  at  1.005  gigahertz 
and  extending  to  8.865  gigahertz  with  a  sampling  interval.  A/,  of  60  megahertz.  These 
frequency  values  provide  a  maximum  unambiguous  range,  Rmaxi  of  2.5  meters  with  a  range 
resolution,  Ar,  of  1.91  centimeters. 


3.8  Results 

3.8.1  Single  Sphere 

The  single  sphere  did  not  provide  a  large  amount  of  information.  While  there  are  a 
number  of  scattering  modes  to  examine,  the  dynamic  range  between  the  specular  return 
and  the  creeping  waves  makes  even  the  first  creeping  wave  difficult  to  see,  and  the  second 
time  around  creeping  wave  is  not  discernable  without  special  processing.  This  second  time 
wave  was  only  visible  using  double  precision  math  for  all  calculation  along  with  increasing 
the  frequency  sampling  interval  and  the  band  of  frequencies  swept. 

Spectrum 

The  initial  spectrum  for  a  sphere  starts  low  and  oscillates  through  the  resonance  region 
as  it  settles  towards  the  final  value  (figure  3.4).  These  oscillations  are  due  to  constructive 
and  destructive  interference  between  the  reflection  from  the  front  of  the  sphere  and  the 
creeping  wave  travelling  around  the  sphere.  The  oscillations  decay  since  the  creeping  wave 
intensity  decreases  with  frequency.  At  first  glance,  the  phase  response  is  essentially  linear 
as  the  number  of  wavelengths  to  the  front  of  the  sphere  increases  with  frequency.  However, 
close  examination  shows  small  deviations  due  to  contributions  from  the  creeping  wave. 

Inverse  Fourier  Transform 

Applying  an  inverse  Foiurier  transform  and  scaling  the  horizontal  axis  by  the  two-way 
propagation  time  provides  a  range  profile  that  clearly  shows  the  leading  edge  of  the  sphere 
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Single  Sphere 


Spectrum 


Figure  3.4:  Single  Sphere:  Spectral  Response 


0.15  meters  in  front  of  the  phase  reference  (see  figure  3.5).  This  plot  also  clearly  shows 
the  creeping  wave  at  its  expected  positions  of  — 7rr/2  =  -0.236  meters.  The  second-time 
around  creeping  wave  should  occupy  a  position  at  — 37rr/2  =  —0.707  meters,  but  it  is  not 
visible  in  this  plot  due  the  extremely  small  contribution  it  makes  over  this  frequency  range. 

Spectrogram 

The  spectrogram  for  the  single  sphere  is  shown  in  figure  3.6.  Here  the  specular  return 
shows  as  a  very  clear  band  located  0.15  meters  ahead  of  the  phase  reference,  while  the 
creeping  wave  is  located  at  an  approximate  position  of  0.24  meters  behind  the  phase  refer¬ 
ence.  The  creeping  wave  component  is  fairly  clear  at  the  lowest  frequencies  and  becomes 
less  pronounced  as  the  frequency  increases.  Once  again,  the  second-time  around  creeping 
wave  lacks  enough  energy  to  be  clearly  discernible,  but  the  distribution  does  exhibit  a  very 
slight  line  at  the  correct  range. 

Wigner  Distribution 

The  Wigner  distribution  for  the  single  sphere  is  shown  in  figure  3.7.  Like  the  spec¬ 
trogram,  the  Wigner  distribution  clearly  shows  the  specular  and  the  creeping  wave  re- 
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Single  Sphere  -  IFT 
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Figure  3.5:  Single  Sphere:  Inverse  Fourier  Transform 


Single  Sphere  Spectrogram 
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Figure  3.6:  Single  Sphere:  Spectrogram 
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Single  Sphere  Wigner 


Frequency  (QHz) 


Figure  3.7;  Single  Sphere:  Wigner  Distribution 


turn.  The  major  difference  is  the  appearance  of  very  large  cross  terms  between  the  ac¬ 
tual  signal  components.  The  actual  signal  component  should  be  the  specular  reflection  at 
-f0.15  meters,  the  creeping  wave  at  —0.235  meters,  and  the  second-time  creeping  wave  at 
—0.707  meters.  These  components  give  rise  to  cross  terms  at  (0.15  —  0.235)/2  =  —0.0425, 
(0.15  —  0.707)/2  =  —0.278,  and  (—0.235  —  0.707)/2  =  —0.471  meters.  This  transform  is 
periodic  along  the  frequency  axis,  so  additional  cross  terms  appear  due  to  the  repetitions  of 
the  three  signal  components.  The  specular  and  the  specular’s  periodic  copy  cause  the  large 
cross  term  which  appears  at  —1.175  meters. 

The  creeping  wave,  although  smaller  than  a  nearby  cross  term,  shows  three  expected 
features.  In  addition  to  appearing  at  the  correct  range,  the  component  has  a  decrease  in 
intensity  as  the  signal  frequency  rises,  and  the  creeping  wave  line  also  has  a  slight  upward 
slope,  indicating  the  propagation  velocity  for  that  mode  varies  with  frequency,  slowing  down 
slightly  as  the  frequency  increases. 

Binomial  Distribution 

The  binomial  distribution,  shown  in  figure  3.8  shares  some  characteristics  with  the 
Wigner,  but  contains  significantly  reduced  cross  terms.  While  the  cross  terms  between  the 
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Single  Sphere  Binomial 


Figure  3.8:  Single  Sphere:  Binomial  Distribution 


specular  and  the  creeping  wave  which  appears  at  —.043  meters  still  exists,  there  are  no 
visible  cross  terms  due  to  the  periodic  images  of  the  components.  This  representation  also 
shows  the  changing  intensity  and  position  of  the  creeping  wave  as  seen  with  the  Wigner. 
Missing  from  this  distribution  is  any  visible  indication  of  the  second-time  creeping  wave. 

Single  Sphere  Summary 

Since  the  variation  in  return  with  respect  to  frequency  was  small,  the  spectrogram 
actually  provided  a  relatively  sharp  view  of  the  target  in  the  joint  frequency-time  domain. 
It  was  unable  to  discern  the  changing  propagation  velocity  in  the  creeping  wave  mode  which 
was  evident  in  the  Wigner  and  binomial  distributions. 

3.8.2  Two  Spheres:  In-Line 

When  positioned  one  behind  the  other,  two  spheres  should  provide  two  specular  scatter¬ 
ing  modes  and  two  readily  discernible  creeping  wave  modes.  These  spheres  will  also  scatter 
based  on  the  interaction  of  the  spheres.  The  primary  interaction  between  the  spheres  is  the 
wave  that  reflects  off  the  back  sphere,  travels  back  to  the  front  sphere  where  it  reflects  off  the 
rear  surface  of  the  sphere,  then  reflects  off  the  back  sphere  a  second  time  before  returning 
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to  the  radar.  Placing  the  phase  reference  at  the  midpoint  between  the  two  spheres  means 
the  specular  returns  should  appear  at  distances  of  +0.65  and  —0.35  meters,  the  creeping 
wave  returns  at  +0.265  and  —0.735  meters,  and  the  triple  reflection  return  at  —1.05  meters. 


Spectrum 


Figure  3.9  gives  the  spectral  return.  The  spectral  response  shows  rapid  variations  in  the 
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Figure  3.9:  In-Line  Spheres:  Spectral  Response 


scattered  signal  as  the  components  combine  in  and  out  of  phase  with  changes  in  wavelength 
of  the  illuminating  plane  wave.  With  the  single  sphere,  these  oscillations  settled  out  as  the 
frequency  increased  because  only  one  significant  mode  existed.  With  the  two  spheres,  the 
creeping  wave  contribution  also  drops  off  with  increasing  frequency  but  the  specular  returns 
continue  to  make  significant  contributions  throughout  the  frequency  sweep. 

Inverse  Fourier  Transform 

The  inverse  Fourier  transform  for  this  signal  shows  five  clearly  discernible  peaks  (fig¬ 
ure  3.10).  The  first  four  peaks  correspond  to  the  speculars  off  the  front  of  each  sphere  (at 
0.65  meters  and  —0.35  meters),  the  creeping  wave  around  each  sphere  (at  .264  meters  and 
—.736  meters).  The  return  at  —1.05  meters  is  due  to  part  of  the  specular  return  off  the 
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In-Line  Sphsree  -  IFFT 


Figure  3.10:  In-Line  Spheres:  Inverse  Fourier  Transform 


second  sphere  being  reflected  by  the  first  sphere  before  scattering  back  to  the  source. 
Spectrogram 

Viewed  with  the  spectrogram  (figure  3.11),  the  two  speculars  and  two  creeping  waves 
are  clearly  visible.  The  doubly  reflected  mode  also  is  visible.  As  the  frequency  increases, 
there  is  a  slight  drop  in  intensity  for  the  creeping  wave  modes.  The  apparent  position  for  the 
creeping  wave  also  moves  slightly  over  the  frequency  range  since  the  creeping  wave  travels 
around  the  spheres  faster  at  lower  frequencies. 

Wigner  Distribution 

The  Wigner  distribution  for  this  signal,  shown  in  figure  3.12,  suffers  from  many  large 
interference  terms.  All  five  components,  two  speculars,  two  creeping  waves,  and  the  triple 
bounce,  along  with  their  periodic  copies  combine  to  generate  cross  terms.  All  this  interfer¬ 
ence  makes  it  difficult  to  directly  view  the  lower  energy  scattering  components. 

Binomial  Distribution 

The  binomial  distribution  shown  in  figure  3.13  illustrates  how  much  clearer  a  reduced 
interference  distribution  can  be.  In  this  case,  most  of  the  cross  terms  that  cluttered  the 
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Figure  3.12:  In-Line  Spheres:  Wigner  Distribution 
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In-Line  Binomial 


Figure  3.13:  In-Line  Spheres:  Binomial  Distribution 


Wigner  distribution  have  been  suppressed.  The  specular  reflections  are  clearly  shown,  and 
the  creeping  waves  off  each  sphere  can  also  be  seen.  The  triple  bounce  component  at  —1.05 
meters  is  too  small  to  show  up  in  this  distribution. 

In-Line  Sphere  Summary 

With  this  geometry,  the  main  contribution  comes  from  the  front  of  each  sphere.  The 
creeping  waves  of  the  two  spheres  provide  a  small  return  best  seen  in  the  spectrogram.  The 
spectrogram  also  shows  the  doubly  reflected  component  at  —1.05  meters. 

3.8.3  Two  Spheres:  Broadside  —  Perpendicular  Polarization 

Positioning  the  two  spheres  broadside  to  the  radar  changes  the  number  and  types  of 
modes  expected  in  the  returned  signal.  In  addition  to  the  returns  from  a  single  sphere, 
several  modes  exist  to  transfer  energy  between  the  spheres  before  sending  it  back  towards 
the  radar. 

The  signal  polarization  also  starts  to  play  a  major  role  in  this  geometry.  Unlike  the 
single  sphere  and  the  in-line  spheres,  this  problem  is  no  longer  axially  symmetric  with 
respect  to  the  direction  of  propagation.  The  first  case  considered  will  orient  the  electric 
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field  perpendicular  to  the  line  connecting  the  sphere  centers.  This  means  creeping  waves 
will  not  scatter  energy  firom  one  sphere  in  the  direction  of  the  other.  All  sphere  interactions 
with  this  polarization  are  based  on  reflections  as  governed  by  geometric  optics.  Each  sphere 
will  act  independently  to  provide  the  creeping  wave  component  as  would  be  seen  for  a  single 
sphere. 

Spectrum 

As  shown  in  figure  3.14,  the  spectral  response  is  similar  in  some  ways  to  the  single  sphere 

Broadside  -  Perpetxiicular  Pol  Spectrum 


Figure  3.14:  Broadside  —  Perpendicular  Polarization:  Spectral  Response 

response.  The  reflections  off  the  front  and  the  creeping  waves  off  each  sphere  reinforce  each 
other  to  provide  a  reinforced  single  sphere  type  return.  The  reflections  travelling  between 
the  spheres  cause  this  spectrum  to  deviate  from  the  single  sphere  response. 

Inverse  Fourier  Transform 

Figure  3.15  shows  the  inverse  Fourier  transform  for  the  broadside  spheres.  The  main 
peaks  correspond  to  the  front  of  the  spheres  at  +0.15  meters,  the  creeping  waves  at  —0.24, 
and  a  reflected  mode  that  boimces  once  off  each  sphere  before  returning  to  the  radar  at 
—0.64  meters. 
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Broadside  -  Perpendicular  Pol  -  I  FT 


Figure  3.15:  Broadside  —  Perpendicular  Polarization:  Inverse  Fourier  Ttansform 


Spectrogram 

The  spectrogram  in  figure  3.16  shows  the  same  three  scattering  modes  discussed  in  the 
previous  section  due  to  the  front  of  the  spheres,  the  creeping  waves,  and  a  double  reflection 
return.  In  addition  to  showing  the  position  for  these  returns,  the  drop  in  creeping  wave 
intensity  as  the  signal  frequency  increases  is  also  evident. 

Wigner  Distribution 

Once  again,  the  Wigner  distribution,  show  in  figure  3.17  shows  the  strong  return  off  the 
front  of  the  spheres  very  clearly,  but  it  has  trouble  with  losing  the  weaker  components  in 
the  cross  terms.  The  creeping  wave  shows  up  at  the  low  end  of  the  frequency  sweep,  but  it 
drops  ofiF  quickly.  The  cross  terms  between  the  front  of  the  sphere  and  the  creeping  wave 
returns  are  evident  across  the  entire  frequency  range. 

Binomial  Distribution 

The  reduced  cross  term  levels  of  the  binomial  distribution,  shown  in  figure  3.18,  make 
it  much  easier  to  view  the  creeping  wave  behavior.  In  particular,  the  drop  in  intensity  and 
the  change  in  apparent  position  are  both  evident  in  the  binomial  distribution. 
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Figure  3.16:  Broadside  —  Perpendicular  Polarization:  Spectrogram 


Broadside  Perpendicular  Wgner 


Figure  3.17:  Broadside  —  Perpendicular  Polarization:  Wigner 
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Broadside  Perpendicular  Binomial 


Figure  3.18:  Broadside  —  Perpendicular  Polarization:  Binomial 


Broadside  With  Perpendicular  Polarization  Summary 

When  illuminated  with  this  polarization,  the  broadside  spheres  acted  very  similarly 
to  a  single  sphere.  Although  modes  existed  which  representated  interactions  between  the 
two  spheres,  their  contribution  was  relatively  small,  especially  when  viewed  in  the  time- 
frequency  domain. 

3.8.4  Two  Spheres:  Broadside  —  Aligned  Polarization 

Keeping  the  same  target  position  and  orientation,  the  radar  system  was  rotated  ninety 
degrees  so  the  incident  electric  field  was  in  alignment  with  the  line  between  the  spheres’ 
centers.  This  polarization  will  support  all  the  scattering  modes  seen  in  previous  section  but 
will  generate  additional  modes  using  creeping  waves  to  transfer  energy  between  the  spheres. 

Spectrum 

The  spectral  response  in  figure  3.19  is  very  similar  for  the  broadside  spheres,  with  the 
phase  response  virtually  identical.  Once  again,  the  magnitude  of  the  spectrum  shares  some 
characteristics  with  the  single  sphere  but  behaves  differently  as  the  frequency  increases. 
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Broadside  -  Aligned  Pol 


Spectrum 


Figure  3.19:  Broadside  —  Aligned  Polarization:  Spectral  Response 


Inverse  Fourier  Transform 

The  inverse  Fourier  transform  show  in  figure  3.20  contains  the  same  peaks  seen  with  the 
other  polarization.  It  shows  the  front  of  the  spheres,  the  creeping  wave,  and  the  reflection 
from  one  sphere,  off  the  second,  and  back  to  the  radar.  Three  additional  modes  appear 
in  this  plot.  The  return  at  —0.53  meters  represents  energy  reflecting  from  the  first  sphere 
toward  the  second,  then  traveling  around  the  back  of  the  second  sphere  before  returning 
to  the  radar.  The  second  new  mode  at  —0.74  meters  follows  a  creeping  wave  90  degrees 
around  to  the  back  of  the  first  sphere,  travels  to  the  second  sphere,  and  reattaches  as  a 
creeping  wave  for  another  90  degrees  before  heading  back  to  the  r^ldar.  The  final  new  mode 
employs  a  creeping  wave  around  the  first  sphere,  then  reflects  off  the  second  sphere  back 
towards  the  first  where  it  reflects  back  towards  the  radar.  This  mode  appears  at  a  range  of 
—0.88  meters. 

Spectrogram 

Figure  3.21  shows  the  spectrogram  for  the  broadside  spheres  with  the  incident  wave 
aligned  with  the  line  between  the  spheres.  Several  components  appear  in  this  distribution 
at  ranges  of  0.15,  —0.24,  —0.35,  —0.53,  and  —0.63  meters. 
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Broadside  -  Aligned  Pol  -  I  FT 


Figure  3.20:  Broadside  —  Aligned  Polarization:  Inverse  Fourier  Transform 


Broadside  Aligned  Spedrogram 


Figure  3.21:  Broadside  —  Aligned  Polarization:  Spectrogram 
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Wigner  Distribution 

The  large  number  of  scattering  modes  leads  to  many  cross  terms  in  the  Wigner  distribu¬ 
tion,  shown  in  figure  3.22  The  only  clearly  illustrated  modes  are  the  specular  reflection  off 
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Figure  3.22:  Broadside  —  Aligned  Polarization:  Wigner 

the  front  of  the  spheres  and  the  creeping  waves.  There  is  a  slight  indication  of  the  multiple 
bounce  reflective  component  at  —0.64  meters. 

Binomial  Distribution 

The  binomial  distribution  in  figure  3.23  does  not  have  large  cross  terms,  but  only  three 
modes  show  clearly:  the  sphere  fronts  at  0.15  meters,  the  creeping  waves,  initially  at  —0.24 
meters,  and  the  combined  specular /creeping  wave  mode  at  —0.53  meters.  The  binomial 
does  clearly  illustrate  the  dispersive  nature  of  the  creeping  wave  mode. 

Broadside  Aligned  Polarization  Summary 

As  expected,  more  interactions  were  evident  between  the  two  spheres  using  this  polar- 
ization.  This  was  due  to  the  direction  of  travel  for  the  creeping  waves  around  the  spheres. 
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Broadside  Aligned  Binomial 
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Figure  3.23:  Broadside  —  Aligned  Polarization:  Binomial 


3.9  Summary 

Applying  Cohen’s  class  distributions  to  generate  time-frequency  distributions  using 
swept  frequency  signals  was  successful  and  informative.  This  technique  provided  informa¬ 
tion  on  the  targets  and  their  scattering  modes  not  readily  seen  using  either  the  frequency 
or  the  time  domain  representations  of  the  signal.  By  using  time-frequency  techniques,  dis¬ 
persive  effects  were  highlighted,  even  though  the  dispersion  involved  with  these  targets  was 
small. 


44 


CHAPTER  4 


STATISTICS  FOR  COHEN’S  CLASS  TFD  OF  WHITE 

GAUSSIAN  NOISE 


4.1  Introduction 

Time-frequency  analysis  is  based  on  transforming  a  signal  from  either  a  function  of 
time,  g{t),  or  frequency,  G{f),  into  a  function  of  both  time  and  frequency,  Cg{t,f).  This 
representation  has  the  advantage  of  capturing  the  features  of  non-stationary  signals  and 
has  been  heavily  used  with  acoustic  signals  [1]. 

Two  frequently  used  time-frequency  transformations  generate  the  Wigner  distribution 
and  the  spectrogram  [37].  These  distributions  belong  to  a  larger  class  of  distributions 
referred  to  as  Cohen’s  class  [12].  The  derivation  here  will  focus  on  this  general  case. 

This  chapter  examines  the  mean  and  variance  for  the  time-frequency  distribution  of 
a  noise  signal.  Both  analytical  and  numerical  values  are  presented  for  a  sampled  signal 
processed  with  discrete  transformations. 

This  analysis  assumes  a  complex  noise  only  signal,  g{t),  comprised  of  two  independent, 
identically  distributed  random  variables,  gr{t)  and  gi{t),  which  are  normally  distributed 


with  a  white  spectra.  This  means  g{t)  is  independent  of  gir)  unless  t  —  t.  The  mean, 
variance,  autocorrelation,  and  spectrum  for  the  noise  signal  are  defined  as: 

mt)]  =  0  (4.1) 

var[5(f)]  =  (4.2) 

Rgi^)  =  (4.3) 

\Gif)f  =  |^b«)]|"  =  a2  (4.4) 

This  white  noise  assumption  creates  problems  with  the  continuous  signal  since  it  implies 


infinite  signal  bandwidth  and  infinite  power.  A  finite  result  is  possible  when  using  a  kernel 
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with  finite  extent  in  the  ambiguity  domain,  but  many  transforms,  such  as  the  Wigner, 
binomial  and  spectrogram,  do  not  satisfy  this  constraint.  However,  the  band-limited  nature 
of  a  discrete  signal  will  allow  this  for  any  transformation  kernel. 

Throughout  this  chapter,  all  integrations  will  have  limits  going  from  —  oo  to  -boo  unless 
otherwise  specified.  With  the  discrete  formulas,  the  indices  run  from  0  to  iV  —  1,  where  N  is 
the  number  of  data  points  considered  or  0  to  P  —  1,  where  P  is  the  length  of  the  transform 
used. 

4.2  Continuous  Transform 

A  continuous  time  signal,  g{t),  can  be  represented  in  the  joint  time-frequency  domain 
using  a  transformation  belonging  to  Cohen’s  class  [12],  expressed  as 

Cg{t,f)  =  1 1  jg(u  +  0  g*  -  0  4>{v,  (4.5) 

The  three  integrations  in  the  preceding  equation  represent  two  forward  and  one  inverse 
Fourier  transforms.  Using  the  relationships 


^tf  [^(^)] 

=  j  x{r)e  dr 

(4.6) 

=  x{f) 

(4.7) 

=  j  X{\)ej^^^^d\ 

(4.8) 

=  x{t) 

(4.9) 

allows  re-writing  the  transformation  as 

Cg{tJ)=P,f 

-u)  5  («  +  0  9*  -  0  j  0(i^,t)0 

(4.10) 

Cohen’s  class  transformations  cover  a  wide  range  of  time-frequency  distributions,  in¬ 
cluding  the  spectrogram,  the  Wigner-Ville,  and  the  reduced  interference  distribution  (RID) 
[29]. 

Prom  the  view  of  implementation,  this  form  involves  three  Fourier  transformations  and 
a  simple  multiply.  Since  multiplies  in  one  domain  correspond  to  convolutions  in  the  other 
domain,  equation  4.10  can  be  re-written  as 

Cgit,  /)  =  /  /  ^  +  0  ”  i)  (4.11) 

or,  in  a  more  compact  form 

Cg{t,  f)  =  Prf  [Rg{t,  r)  0t  tpit,  r)]  (4.12) 
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with  Rg{t,T)  as  the  local  autocorrelation,  representing  the  convolution  with  respect  to 
t,  and  as  the  autocorrelation  domain  kernel. 


4.2.1  Continuous  Mean 

Performing  a  Cohen’s  class  transformation  on  this  signal  gives  a  time-frequency  distri¬ 
bution  Cg{t,f).  Start  with  the  expected  value  of  this  distribution 

y  j  j  3  (^  +  ^9*  (u-  0  (4.13) 

where  the  operator  £^(t,/) [Cg(t,  /)]  represents  the  expectation  of  Cg{t,  f)  and  can  be  written 
explicitly  as 

E(t,f)  [Cg{t,  /)]  =  //  Cgit,  f)fg{t,  f)dtdf  (4.14) 

using  fg{tyf)  to  represent  the  probability  distribution  function  for  the  noise  signal. 

The  expectation  operation  can  move  inside  the  integrals  over  r  and  u  since  they  represent 
linear  transformations  between  time  and  frequency  domains.  This  gives 


E^tJ)[Cgit,f)]  = 

/  /  E(,,,r)  [J  9  +  0  9*  0  -  0  e^^^‘''^du<j){u,T)  (4.15) 

The  expectation  operation  will  return  zero  unless  r  =  0.  For  the  case  of  t  =  0,  the 
expectation  is  a^.  This  allows  rewriting  the  equation  as 

E[Cg{t,f)]  =  I  j  j  oH{T)4>{v,ry‘^^^'''^-'^-f'^UududT  (4.16) 

Moving  the  constant  outside  the  integral  and  nesting  the  three  integrals  yields 

(t,  /)]  =  (/  6{T)e-^^^f^dT  (4.17) 

The  inner  integration  gives  the  Fomier  transform  of  a  constant  which  leads  to  a  delta 
function  in  the  frequency  domain,  S{u). 

E[Cg{t,f)]=a^  j  J  6{i')(p{v,T)e~^‘^^''*di'6{T)e~^^^^'^dT  (4.18) 

Both  integrals  contain  a  delta  function  of  the  integration  variable.  This  means  the 
integrand  is  sampled  at  the  point  where  the  argument  of  the  delta  function  equals  zero. 


Cgit,f)=E[Cg{t,f)]  =  a^4>{0,0) 


(4.19) 
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The  expected  value  is  the  noise  power  (the  variance  of  the  noise  signal,  or  a^)  multiplied 
by  the  transformation  kernel,  evaluated  at  the  origin.  For  the  distributions  of 

concern,  (j)  has  a  value  of  1  at  the  origin  [7],  giving  the  final  result 

E[Cg{t,f)]  =  a^  (4.20) 

This  can  also  be  evaluated  using  the  autocorrelation  expression  of  the  transform,  which 
begins 

E[Cg{t,  f)]  =  j  j  E^g(u  +  0  g*  (u  -  0  ip{t  -  u,  (4.21) 

Once  again,  the  expectation  returns  a  non-zero  value  only  if  r  =  0  in  which  case  it  becomes 
a^.  Substituting  this  relationship  and  reversing  the  order  of  integration  means 

E[Cg{t,f)]  =  j  J  a^6{T)‘ip{i  —  u,T)e~^^'^^'^dudT  (4.22) 

=  J  i]}{t  —  u,  0)du  (4.23) 

=  (4.24) 

The  last  step  in  the  two  preceding  derivations  assume  the  kernel  generates  a  distribution 
which  satisfies  the  time  marginal  property.  To  satisfy  the  time-marginal  property,  </>(/,  0) 
must  be  one  [22],  which  means  the  kernel  in  the  autocorrelation  domain  must  have  the 
property  ip{t,  0)  =  6{t). 

4.2.2  Continuous  Variance 

Finding  the  variance  for  the  time-frequency  distribution  is  more  complex.  The  derivation 
begins  with 

var  [Cgit,  /)]  =  E  [\Cg{t,  f)  -  E[Cgit,  /)]|2]  (4.25) 

Since  the  expected  value  of  Cg{t,  f)  is  a  constant,  the  argument  of  the  expectation  can 
be  expanded  and  simplified. 

var[Cg(t,/)]  =  E[Cg{t,f)C;{t,f)-Cg{t,f)C;^ 

-  c;{t,  f)C^ -f-  q(U)  ]  (4.26) 

=  E[Cg{tJ)C;{tJ)]-E[Cg{tJ)]C^ 

-E  [c;  {t,  /)]  cjifj) 

=  E[Cg{t,f)c;it,f)] -Cg{t,f)  c;it,f)^ 

=  E[Cgit,f)C*git,f)]-a^mO)\^ 
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(4.27) 

(4.28) 

(4.29) 


The  first  term  can  be  expanded  using  the  definition  for  the  transformation 


C,(t,f)  =  J  j  Jg(u+  0  g'  (u  -  0  4,{„,  (4.30) 

which  means 


b[|c,((,/)P]  = 

E\j  //s  0  5*  “  0 

■///«■  (”+ 5) «(”  -  0  (f>*{X,s)e>^^^-^^+^*+^^UvdXdsj 


(4.31) 


Reordering  the  integrals  and  moving  the  expectation  to  include  only  the  nondetermin- 
istic  elements 


E 


\Cg{t,f)\' 


//////®  K“+ i)  «■  (“  -  i) + 1) « (”  - 1). 

T)cf>*  (A,  (4.32) 


Evaluate  this  integral  using  the  expectation  properties  of  white  gaussian  noise.  The 
values  of  g{ti)  and  g{t2)  are  independent  unless  ti  =  *2-  Independence  allows  separating 
the  expectation  of  their  product  into  the  product  of  their  individual  expectations.  For  the 
case  ti  t2 

E[g{ti)g{t2)]  =  E\gih)]E\g{t2)]  (4.33) 

When  applied  to  equation  4.32,  if  ti  ^  t2  ^  h  ^  t4,  then  the  expectation  is 

E[giti)g*it2)g*ih)gih)]  =  E[gih)]E\g{t2)]E\g{t3)]E\g{U)]  (4.34) 

=  0  •  0  •  0  •  0  (4.35) 

=  0  (4.36) 

since  the  individual  distributions  are  zero-mean  gaussian  random  variables. 

A  non-zero  value  results  unless  all  four  times  are  matched  in  pairs,  that  is,  only  two 


time  instances  ti  and  t2  (which  may  be  the  same),  need  to  be  considered 

■E^b(^i)5*(ii)5*(*i)5(ii)]  =  •£' [l5(ii)l^]  (4-37) 

Blg(k)g-{h)g-{h)g{t2)]  =  b[|9(<i)P]b[|9(*2)P]  (4.38) 

£[9(*l)9‘(«2)9‘(*l)s(f2)l  =  £  [!<;((, )P]b[|9(*2)|"]  (4.39) 

B[9(tl)9*((2)!)*(i2)9(tl))  =  B[{9(*.)}"]£:[{!l'(f2)}"]  (4.40) 
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Equation  4.37  requires  u  —  v  with  the  additional  condition  that  t  and  s  are  both  zero. 
Expanding  this  for  the  complex  value  Qri  +  jgn  gives 


E[\9rl  +^1^] 

(4.41) 

E  ffrl  +  511  +  ^9rl9il 

(4.42) 

3a f  -1-  3cr^  -|-  2a^af 
a^  a^  ^a^  a"^ 

(4.43) 

(4.44) 

2ct^ 

(4.45) 

since  gr  and  gi  have  the  same  distribution  the  total  noise  power  is  split  evenly  between  the 
real  (in-phase)  and  imaginary  (quadrature)  components,  making  —  (t^/2  and  =  (7^/2. 

Equation  4.38  corresponds  to  u  v  and  r  =  s  =  0.  Expand  using  gri  +  jgn  and 

gr2  +  jgi2 

r  r^"ir  r  ^-11- 

(4.46) 

(4.47) 

(4.48) 

(4.49) 


E[\g{t,)f]E[\g{t2)f 

—  E  +  jgn  1^ 

E  |^|5r2  +j9i2f 

=  E  [sri  +  5ii]  E 

9h+9i2 

=  +  (Tiffi  + 


=  a 


r^T 
A 


The  same  result  occurs  for  equation  4.39,  used  when  u  =  v  and  r  =  s  0. 

Equation  4.40  is  valid  when  u  =  v  and  r  =  -s  0.  Here  the  values  at  each  time 
instance  do  not  get  conjugated  and  multiplied  to  give  the  magnitude  squared. 


E  [{5(^1)}^]  E  [{^*(^2)}^  =  E  [(3rl  +jgil)^  E  [(gr2  -jgi2)^]  (4.50) 

=  E  ~  9il  +  2j5rl5il]  E  ^g^2  ~  9i2  ~  '^j9T29i2^  (4-51) 
=  -  (^i  +  o)  ((^r  -  (^i  +  o)  (4.52) 


_  _4  I  4  cy_2_2 

““  (Tj.  I  ^ ^(7 ^ % 

=  0 

with  Or  =  =  ol\/2  allowing  the  move  from  equation  4.53  to  equation  4.54. 

The  three  non-zero  cases  can  be  accounted  for  by  the  substitution 


(4.53) 

(4.54) 


^  K“  i) ‘  0  (”+ f  )"('’■  0. 


=  [5(t)^(s)  -I-  6{u  —  v)6(t  —  s)] 

(4.55) 


into  equation  4.32 
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JUJU  [6{r)6{s)  +  6{u  —  v)S{t  —  s)] 

■(t>{u,  T)(t>*  (A,  s)e^^^^‘^-’'^-f'^-^'’+^*-^^‘‘'>dudi^dTdvdXds  (4.56) 


Rewrite  this  integral  as  a  sum  of  two  terms,  one  accounting  for  the  case  t  =  s  =  0  and 
the  other  r  =  s  ^  0. 

E[\Cg{t,fU]=a\A  +  B)  (4.57) 

where 


mil U(t^*  s)e^^^^''^-‘'^-f^-^^+^^+f^UududTdvdXds  (4.58) 


-  =  /////^-  v)6{t  -  s)(l)iu,  t)(P*{X,  s)e^^^^‘''^-‘'^-f^-^^+^^+f‘>UududTdvdXds 


(4.59) 


In  A,  the  expression  6{s)6{r)  allows  removing  two  levels  of  integration  resulting  in 

A  =  Jjl  j<f>{u,0)<l)*{X,0)e^^^<^‘'^-‘'*-^'’+^*UududvdX  (4.60) 

Regrouping  to  isolate  the  u  and  v  dependencies 

0)e^^^^^^-‘^UiydX  (4.61) 

The  inner  integrals  act  as  Fourier  transforms  of  constant  arguments  and  can  be  replaced 
with  the  delta  functions,  S{u)  and  5(A),  leaving 

A  =  j  j  6iu)6{X)(l>{u,0)<l>*  {X,0)ej^^^^^-'^UudX  (4.62) 

which  simplifies  to 

A-0(O,O),^*(O,O)  =  |</.(O,O)|2  (4.63) 

In  the  expression  for  B,  6{u  —  v)  and  5(r  —  s)  allow  the  substitution  of  v  =  «  and  s  =  r 
and  removal  of  the  v  and  s  integrals,  leaving 


//(/ 


^P^'''^du 


B  =  jjj  j  (l>{u,T)<l>*{X,T)e>^<-^^+‘^+^^-^^UududTdX  (4.64) 


Evaluating  the  integration  over  u  results  in 


=  5(i/  —  A) 


(4.65) 

(4.66) 
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When  substituted  into  the  expression  for  B,  this  allows  removing  the  A  integration.  B 
becomes 

B  —  j  j  (4.67) 

~  1 1  (4.68) 

Taking  the  expression  for  variance  in  equation  4.29  and  substituting  equations  4.63  and 

4.68  {A  and  B) 

var[C'g(t,/)]  =  E[\Cg{t,f)f]-E[Cg{t,f)f  (4.69) 

=  +  -CT^|(/»(0,0)p  (4.70) 

=  a^|(^(0, 0)1^  Tcr^yy' 10(1/,  r)|^dz/dr-(7^|^(0,0)p  (4.71) 

=  J  J  |0(i^)7')Pdi/dr  (4.72) 

Unfortunately,  this  integral  cannot  be  evaluated  for  all  possible  kernels.  A  popular  class 
of  kernel  functions  are  product  kernels  —  that  is,  0(i/,  r)  =  0(i/t)  —  with  the  value  of  1 
at  the  origin  (necessary  for  desirable  marginal  properties  in  the  distribution  [21].  For  these 
kernels,  the  integrals  are  unbounded. 

Viewed  from  the  autocorrelation  domain,  the  mean  square  value  is  represented  by 

E[\Cg{t,f)f]  =  e'^I  lg(u  +  ^'^g*(u-^^^l;{t-u,T)e-^^”f-dudr 

•  y'  fir  +  0  5*  -  0  0(t  -  V,  s}e~^^^f^dvds  (4.73) 

Moving  the  expectation  to  cover  only  the  random  elements  and  reordering  the  integrals 

=  ////^K«+0i>-(-i)s(''+0»*(»-|); 

•0(t  —  n,  T)0(t  —  V,  dudrdvds  (4.74) 

As  described  earlier,  the  expectation  can  be  represented  inside  this  integral  as  two  delta 
fvmctions. 

E[\Cg{t,f)f]  =  j  1 1  l[6{T)6{s)+6iu-v)6{T-s)]<r^ 

•ip^t  —  u,  T)'ip{t  —  V,  dudrdvds  (4.75) 

Writing  this  as  the  sum  of  two  integrals  and  applying  the  properties  of  the  delta  functions 
gives  the  following  result 

E[\Cg{tJ)f]  =  —  u,  0)0 (i  —  V,  0)dudv  +  J  J  ip{t  —  u,  —  m,  —T)dudT 

(4.76) 
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To  simplify  this,  the  standard  kernel  property  =  S{t)  [22]  allows  removal  of  the 

first  integral.  Using  a  change  of  variables,  integration  limits  of  infinity,  and  symmetry 
properties  of  standard  kernels  allows  writing  the  t  —  u  arguments  as  u. 

=  J  'ip{u,T)'il){u,-T)dudT  (4.77) 

The  variance,  obtained  from  subtracting  the  square  of  equation  4.24  from  equation  4.77, 
has  the  same  form  as  found  earlier  when  working  in  the  ambiguity  domain. 

va.v[Cg{t,f)]  —  a^J  j 'tl^{u,r)ip{u,—T)dudT  (4.78) 

As  seen  before,  these  integrals  with  infinite  limits  are  unbounded. 

4.3  Discretized  Transform 

Results  form  the  continuous  transformation  shows  an  ill-behaved  expression  for  variance 
because  of  the  infinite  power  required  for  a  truly  white  noise  signal.  For  discrete  signals,  a 
white  signal  is  physically  realizable  because  discrete  signals  are  inherently  band-limited. 
Returning  to  the  continuous  definition  of  the  time-frequency  distribution 

Cg{t,f)  =  J  J  j  g  (u  +  0  g*  -  0  <j){v,  T)e^‘^'^^''^~''^~'^'"'>dudvdT  (4.79) 

To  discretize  this,  nAT  will  represent  time  andpA/  the  frequency  where  AT  is  the  sampling 
interval  and  A/  the  corresponding  frequency  resolution. 

Cg{nAT,pAf)  = 

mAT  +  g*  (mAT  -  <f>{qAf,  lAT) 

I  q  m  ^  ' 

,^2Tr{gA  fmAT-gAf  nAT-pAflAT)  j  (4.80) 

The  sampling  interval,  AT,  can  be  set  to  1  second  without  loss  of  generality.  Since 
the  frequency  resolution.  A/,  is  given  by  the  reciprocal  of  the  window  length,  NAT,  the 
summations  can  be  re-written  as 

=  ;^  E  E  Es  (™  +  5)  S*  (”•  -  0  (f . ')  (4.81) 

4.3.1  Discretized  Mean 

To  find  the  expected  value  for  the  discrete  distribution,  the  non-deterministic  g{n)  terms 
need  to  be  kept  inside  the  £?[•]  operator 

E[Cg{n,p)]  = 
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^  E  E  E ®  [«  (-  +  0  5*  (”  -  0]  ^  (4.82) 


The  expectation  for  the  gg*  term  is  a^Si,  so  the  equation  becomes 
E lC,(r>,p)]  =  E E E ‘’"‘‘•I' 

I  q  m  '  ' 

Resolving  the  I  summation 

£:[Cj(n,p)l  =  4EE<»(^.o)'^^‘”’ 

q  m  ^  ' 

O'  '  m 


(4.83) 

(4.84) 

(4.85) 

q  N-  '  /  j7i 

The  second  summation,  which  involves  the  complex  exponential  with  respect  to  m  will 
evaluate  to  zero  miless  g  =  0,  in  which  case  the  summation  becomes  iV.  The  mean  value 
may  be  re-written  as 


E  [Cg{n,v)]  =  a2  X:  o)  (4.86) 

Performing  the  last  siunmation  and  applying  the  property  <^(0, 0)  =  1  to  obtain 

E[Cg{n,p)\  =  a^  (4.87) 

which  matches  the  expected  value  found  for  the  continuous  case. 

4.3.2  Discretized  Variance 

Turning  to  the  variance  calculation,  we  begin  with  the  discrete  Cohen’s  class  transfor¬ 
mation  from  equation  4.81.  As  previously  explained,  the  new  term  needed  is 


[|C's(">p)P_  = 


^EEEH»  +  0/(-n-0<^,<)e’^<~'} 

~EEE  {.  (=+  I)  (-  0  0]  (4.88) 

Rearranging  the  summation  symbols  and  collecting  random  variables  under  the  expec¬ 
tation  operator 

E[\Cg{n,p)\^]  = 

i.EEEEEE^K"..0.-(™-0.-(c.0«(c-0] 

•</.  (J^,  (I)*  e^^i<im-bc)+{b-q)n+p{d-l)  (4 
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As  with  the  continuous  variance,  the  expectation  operation  has  three  cases  with  non-zero 
results.  These  are  accounted  for  by  substituting 

E  g{m+^^g*  {m-  0  g*  ^  ^  +  ^{m-c)\l-d)]  (4.90) 

The  variance  becomes 

E[\Cgin,p)\'^^  = 


jY2  ^  ^  ^  +  ^(m-c)'^(J-d))] 


I  q  m  d  b  c 


.(f,  <t>*  eJ^(9m-bc+(6-9)r*+p(d-0) 


(4.91) 


As  with  the  continuous  case,  breaking  this  expression  into  two  pieces  and  incorporating 
each  8  function  relationship  results  in 

E[\Cg{n,p)\^]=^[A  +  B]  (4.92) 


where 


A  = 


I  q  m  d  b  c  \-^'i  /  XJ’  / 

.g;i^{qm-bc+{b-q)n-]^p{d-l)) 


B  = 


I  q  m  d  b  c  /  \^'i  , 


.J^(<im-bc+{b-q)n+p(d-l)) 


(4.93) 


(4.94) 


The  expression  for  A  simplifies  when  summed  over  I  and  d  because  of  the  properties  of 
the  6  functions 


^  =  EEEE  (W’ (4.95) 

q  m  })  c  ^  ^  ' 

E  <!>  o)  f  o)  (4.96) 


=  EEE 

q  m  b  I  c 

The  bracketed  summation  has  two  possible  values.  If  6  /  0,  the  summation  goes  to 
zero.  When  6  =  0,  it  equals  N.  In  that  case 


a  =  nY^ 


(4.97) 


q  im 

The  summation  over  m  will  also  become  N  when  q  =  0  and  is  zero  for  all  other  q.  This 
simplifies  the  expression  for  A  to 


A  =  N^^iO,  0)<j)*  (0, 0)  =  iV2|0(0, 0)p 


(4.98) 
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Applying  the  6  functions  to  simplify  B 


B  =  EEEE  (P  f±j'\  ^j^iQm-bm+bn-gn)  (4  99^ 

I  q  m  b  / 

=  E  E  E  [e  " '’‘H  Hw’  (i- ')  ^ 

I  q  b  I  ^  ' 

The  bracketed  summation  is  non-zero  only  when  q  =  b,  when  it  evaluates  to  N.  This 
leaves  the  result 

B  =  (4-101) 

I  q 

Substituting  for  A  and  B  in  equation  4.92  and  using  the  result  in  the  expression  for 
variance 

var{C3(n,p)}  =  E\^Cg{n,p)\^^- Cg{n,pf  (4.102) 

=  ^[A-fB]-(aV(0,0))2  (4.103) 

I  I  q 

I  q  I 

These  results  parallel  the  continuous  expressions  for  A  and  B  obtained  in  the  preceding 
section. 

In  contrast  to  the  infinite  integral  limits  used  for  the  continuous  case,  the  discrete 
variance  calculation  involves  finite  summations.  This  allows  numerically  evaluating  the 
variance. 

4.4  Sampled  Transform 

A  parallel  approach  considers  evenly  spaced  gaussian  noise  samples,  g{n),  where  g(n)  = 
g{nAT).  Transform  the  g{n)  samples  using  the  discrete  counterparts  to  the  Fourier  trans¬ 
forms 


iV-l 

:Fnp[x{n)]  =  ^  a:(n)e-J^^’"  (4.106) 

n=0 

=  X(p)  (4.107) 

^pniXijp)]  =  (4.108) 

p=0 

=  x{n)  (4.109) 


(4.104) 
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These  allow  expressing  the  transformation  of  the  time  signal  to  the  sampled  time-frequency 
domain  as 


Cg{n,p)  =  Tip  9  ^  g*  (m  -  0  <^(g,  /)|  j  (4.110) 


The  Fourier  transforms  in  equations  4.106  and  4.108  involve  a  periodicity,  either  real  or 
implied,  over  a  period  of  NAt  in  the  time  domain  which  provides  a  periodic  spectra  with 
frequency  resolution  of  Af  =  l/(iVAt),  with  repeats  every  1/At  hertz. 

To  express  the  transform  using  the  autocorrelation  kernel,  convolve  the  kernel  with  the 
autocorrelation  and  take  the  Fourier  transform  with  respect  to  the  lag  variable,  1. 


Cg{n,p)  —  Tip  [Rg{n,l)  (8»„  tp{n,l)]  (4-111) 

where  Rg{n,l)  in  the  local  autocorrelation  and  is  the  autocorrelation  kernel.  The 

Fourier  transform  of  the  autocorrelation  kernel  gives  the  ambiguity  domain  kernel, 

TnqbPin,  /)]  =  <l>{q,  1)  (4.112) 


4.4.1  Data  Limitations 

There  is  an  important  limitation  to  keep  in  mind  when  using  this  sampled  data  approach. 
The  continuous  definition  for  Cohen’s  class  implies  time  shifts  of  ±1/2  when  determining 
the  discrete  autocorrelation  function  Rg{n,l).  For  odd  values  of  I  this  requires  values  at 
points  between  the  samples.  Two  techniques  for  dealing  with  non-integer  sample  points 
yield  quite  different  results  [21]. 

One  approach  treats  these  intermediate  values  as  zero.  Unfortunately  this  traditional 
approach  discards  information  provided  by  odd  lag  values  and  halves  the  maximum  una¬ 
liased  frequency. 

An  alternative  approach  uses  separate  methods  for  odd  and  even  lag  values.  By  slightly 
modifying  the  definition  for  the  kernel  and  the  autocorrelation  function,  it  is  possible  to 
avoid  undersampling  the  signal.  This  alias-free  formulation  allows  an  unambiguous  fre¬ 
quency  range  of  half  the  sampling  frequency  in  comparison  with  the  traditional  formulation, 
which  has  a  range  of  one-fourth  the  sampling  frequency. 

The  derivations  below  follow  the  alias-free  formulation  in  both  ambiguity  and  autocor¬ 
relation  domains.  The  traditional  methodology  differs  by  eliminating  terms  with  odd  lag 
values. 
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4.4.2  Sampled  Mean 


Calculating  the  expected  value  for  the  time-frequency  distribution  computed  from  sam¬ 
pled  data  requires  evaluating 


E[C,{n,p)]  =  E  [/ip  [j  {m  +  0  s’  (™  -  5)]  «9,o}) 


(4113) 


When  expanded,  this  becomes 

P  -t  N  N 


E[Cgin,p)]=E 


1=0^^  q=0m=0  ^  ^ 


g“  J  g- j  ^  qn  ^  qm 

(4: 


14) 


Moving  the  deterministic  portions  of  this  expression  outside  the  expectation  operator 
P  N  N 


1 


/=0  ^=:0m=0 


i 

9  [m  +  -]  9 


("•-0. 


(4.115) 

The  expectation  will  be  non-zero  only  when  1  =  0  and  the  arguments  for  g  and  g*  are 
equivalent.  When  /  7^  0,  the  independence  of  the  samples  drives  the  expectation  to  zero. 
Accounting  for  this  by  replacing  the  E[gg*]  term  with  <7^61,  makes  the  expected  value 

P  N  N 


E[Cgin,p)]  =  ^EE  E 

/=0  gF=0m=0 

Using  the  delta  function  to  resolve  the  summation  over  I 


(4.116) 


E[Cg{n,p)]  =  ^E  E  0(9,O)e-^'^''"e^'^''- 


(4.117) 


^—0  m=0 

In  the  summation  over  m  make  the  substitution 

N 

iZ2L 


J:  eiw”-  =  N6, 


(4.118) 


m=0 


The  result  is 


2  Jv 


ElCg(n,p)]  = 


9=0 


crV(0,0) 


(4.119) 

(4.120) 


To  satisfy  the  marginal  properties,  the  kernel  equals  one  at  the  origin  [18],  with  the  final 
result 

ElCg(n,p)]=a^  (4.121) 
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4.4.3  Sampled  Variance 


Finding  the  variance  of  a  distribution  requires  both  the  mean  value  and  the  mean  square 
value.  The  mean  square  value  calculated  using  the  ambiguity  domain  form  of  the  time- 
frequency  distribution  given  in  equation  4.110  is 

E  [\C,(n,p)f^ 

=  [i,  (m  +  0  g-  (m  -  0]  .f(5,o}) 

■  (^i.  [5  (c+  0  9*  (c-  I)]  «(6,d)})]‘}  (4.122) 


^Cb  9  c-f-U  c--  U' 


Expanding  the  Fourier  and  inverse  Fourier  transforms 


(4.123) 


E[\Cgin,p)\^] 

=  ^  9  (m  +  ^)  g*  (m  - 

1=0-^^  q=Qm=0  ^  ^ 

i:  V  S  ('  +  f)  «  -  5)  (4.124) 

d=0'^^6=0c=0  \  ^/  \  ^/ 

Moving  constants  outside  the  summations  and  deterministic  portions  outside  the  expec¬ 
tation 


E[\Cg{n,p)\^] 

.  P  N  N  P  N  N  r  .  ,v  , 

=  ]^EEEEEE^;9(-"  +  5)9-(m-09-(c+f)9(c-0 

;=0^=0m=0d=0  6=0c=0  ^  \  \  LJ  \  LJ  \  IJ 


)  q=D  m=D  d=[} 

.a./,  -V  «27r. 


■(t>{q,  l)(l>*{b,  d)e^’-FP(d-0ei#(<7-i>)nei7r«’‘^e-if  9'" 


(4.125) 


To  avoid  an  expectation  of  zero,  the  four  samples  under  consideration  must  be  matched 
in  conjugate  pairs.  A  non-zero  value  occurs  when  m  =  c  and  I  =  d,  or  when  I  =  d  =  0. 


m  =  c 

II 

a, 

II 

o 

2a^ 

m  =  c 

m  ^  c 

o 

II 

II 

The  following  substitution  achieves  this  result: 

E  9(m  +  ^g*  (m-  0  g*  ^  ^  (4.126) 
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The  mean  square  value  becomes 


E 


w 


‘  EEEEEEW' 

1=0  q—0  m=0  d—0  6=0  c=0 

+  E  E  E  E  E  E  ^i^d4>{q,  iWib,  )  (4.127) 

;=0  q—O  m=0  rfrrO  6=0  c=0  J 

Resolve  the  summation  over  the  delta  functions 
B[|C,(n,p)|"]  = 

4  (eE  E  E*(9,i)0-(6.i)e'»'«-‘>Vlf(‘-«>’" 

\;=0  g=0m=0  6=0 


g=0m=0  6=0c=0 


(4.128) 


In  the  first  set  of  summations,  the  only  m  dependence  occmrs  in  the  complex  exponential 
term.  This  term  becomes  N  when  q  =  b  and  is  zero  for  all  other  cases. 


E[\Cg{n,p)f  = 


P  N 


+  E  E  EE'^(«>0)<^*(6,0)e^'^(«-*’)”e^  N*’"e-^  N'?^ 

q=0  m=0  6=0  c=0 


(4.129) 


For  the  second  set  of  summations,  c  and  m  dependencies  occur  only  in  the  exponential 
terms  and  have  a  non-zero  result  only  when  9  =  0  and  6  =  0. 

^4  /  P  N 

JV2 


E 

This  reduces  to 


<=0  q=0 


(4.130) 


E[|C'3(n,p)|2]  =^f;f;|<^(g,/)|2+a^  (4.131) 

1=0 q=0 

The  variance  is  the  mean  squared  value  minus  the  square  of  the  mean  value  found  in 
equation  4T21 


var  [Cg{t,  /)]  =  ^  E  E 


(4.132) 


1=0  q=0 
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Determining  variance  in  the  autocorrelation  domain  requires  evaluating 

vax[Cg{n,p)]  =  E  [|Cp(n,p)p]  -  {E[Cg{n,p)\f  (4.133) 

Evaluate  this  expression  by  separating  the  odd  and  even  lag  values,  /,  to  allow  an  alias- 
free  formulation  for  Cg(ra,p)[21]. 


Y  9  (k  +  9*  (k  -  ^{n  -  k,l)e  ^ 


P-1  CX) 

I 

k=-oo 

P-1  CX) 

I 

k=—oo 


-j^pl 


1=0 
even  I 


1=0 
odd  I 


(4,134) 


where  ^'(n, /)  in  the  modified  ambiguity  free-kernel  for  odd  values  of  L 

Modifying  the  lag  index  covers  positive  and  negative  values,  ranging  from  —P/2  +  1  to 
P/2  —  1.  The  local  autocorrelation  for  a  lag  of  —P/2  must  be  zero  to  satisfy  symmetry 
properties  of  the  transform. 

p_i 

Cg{n,p)  =  Y  Y  9  (k  +  ^)  9*  (k  -  ip{n  -  k,l)e-^^P^ 

1=^+1  fc=-oo  \  ^/  \  ^/ 

even  I 


1=^+1  fc=-oo 
odd  I 


(4.135) 


This  expectation  becomes 

E[\Cgin,p)\'^]  =  E 


Y  Y  9(k+^^9*  {k-^^il;{n-kj) 


l==^+\  *=-°o 

even  I 


^-1 


+  Y  Y  9(k  +  ]^+^^9*  (k  +  ]^-^^'<p'{n-k,l)e 


l=^+\  k=-oo 

odd  I 


2t 


(4.136) 

Breaking  up  this  incredibly  cumbersome  mess  into  only  slightly  more  workable  night¬ 


mares 


E[\Cg{n,p)\‘^]  =  E[A  +  B  +  C  +  D] 


(4.137) 
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=  E[A]^-E[B]-\-E[C]-\-E[D] 


(4.138) 


where 


A  = 


^-1 

2  ^  OO  OO 

E 

E  E  E  i 

<-^+1 

k=~ooTn=—oo 

even  1 

even  d 

9  ^m  + 

0  9*  ("i  -  0  -  k 

£—1 

2  ^  00  C30 


B  = 


E  E  E  + 

l^^+l  d=^+l  fc=-oom=-oo  \  V 

even  I  odd  d 

5  ("i  +  ^  +  0  /  ("i  +  ^  -  0  ^(n  -  k,  l)ij}'{n  -  m,  d)e~^ 


--1  ^-1 

2  ^  2  ^ 


2  + 


=  E  E  E 

<=^+1  d==f+l  fc=-oom=-oo  V  ^  IJ  \ 
odd  I  even  d 

■g  +  0  g*  -  0  tp'{n  -  k,  l)rp{n  —  m, 


E  E  E  E  +  ?  +  + 

-P  ,  ^  .  -P  .  .  t— =  V  ^  \  ^ 


D  = 


1=^+1  d=^+l  k=-oom=-oo 

odd  I  odd  d 


•5  +  ^  +  0  5*  +  i  -  0  tp'in  -  k,  l)7p'{n  -  m,  d)t 

The  expected  value  of  A  simplifies  using  the  relationship 


£^[5(*=  +  5)s*(k-5)9-(m+0s(m-0 


=  a'^{Si6d  +  Sk-m^l- 


so 


=  a" 

E 

E 

i=^+i 

-h 

II 

even  1 

even  d 

•ip{n  —  k,  —  m,  ( 

OO 

OO 

= 

E 

E  V’(« 

_k—-oorr, 

1=— OO 

OO 

+ 

E 

E 

1=^+1 

k=—oo 

OO  OO 


even  I 


(4.139) 

^p(i+d)  (4140) 

i) 

(4.141) 

i) 

^p{i+d)  (4142) 

-d)  (4.143) 

] 

(4.144) 

(4.145) 
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Following  the  same  basic  steps  used  to  find  E[A\  will  lead  to  the  expected  value  of  D. 
Since  I  and  d  must  be  odd  in  the  expression  for  D,  they  cannot  be  zero.  This  means  there 
will  not  be  the  6idd  term  as  in  equation  4.144.  The  expected  value  of  D  simplifies  to 

f-l  oo 

E[D\  =  a^  ^  ^  —  kiVjil)' {n  —  k,—l)  (4.146) 

fe=-oo 

odd  I 

The  expected  values  of  B  and  C  are  zero.  In  equation  4.140  the  g{m  +  .5  +  d/2)  and 
g*{m  +  .5  —  d/2)  terms  can  combine  as  a  complex  conjugate  pair  only  if  d  =  0,  which  is  not 
possible  since  the  summation  limits  d  to  odd  values.  The  g{m+.5+d/2)  and  g*{m+.5—d/2) 
cannot  match  up  with  the  g{k  +  l/2)  and  g*{k  — 1/2)  terms.  The  spans,  |d|  and  |/|,  cannot  be 
equal  since  I  must  be  even  and  d  must  be  odd.  The  same  limits  on  the  summation  variables 
force  the  value  of  C  to  zero  in  equation  4.141. 

The  variance  of  the  distribution  now  becomes 


var[C'3(n,p)]  =  E[\Cgin,p)\‘^]- \E[Cg{n,p)]\‘^ 

=  E[A]  +  E[D]-iay 

I  OO  OO 

Z]  H  ip{n-k,0)i){n-m,0) 


(4.147) 

(4.148) 


=  a 


k=—00  77l=  — oo 
^-1 

o  no 


+  Z]  Z  fp{n-k,l)rp{n-k,-l) 

even  I 


1=^+1  k=-oo 


“T-l 


+  Z  Z  ‘>P'{n-k,l)ip'{n-k,-l)-l 

1=^+1  k=-oo 


odd  I 


(4.149) 


Since  all  the  k  based  summations  have  infinite  limits,  the  specific  value  of  n  does  not 
matter,  allowing  n  =  0. 

vax[Cg{n,p)] 


oo  oo 


f-l 


■1+  Z  Z  ^{-k,0)iP{-m,0)  +  ^{-kMi-k,-l) 

kz=:  —  00  m=  — OO 

even  I 


i==f+i  *=-°° 


p 
2  ^ 


+  Z  Z  i’'{-k,m'{-k,-i) 

odd  I 


l==f+l 


(4.150) 
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If  we  assume  the  kernel  provides  a  time-frequency  distribution  with  the  correct  time 
marginal,  the  kernel  at  zero  lag  can  be  treated  as  a  delta  function,  ^(n,  0)  =  Sn,  forcing  the 
first  summation  to  a  value  of  1 


var[C'g(n,p)] 


=  a 


^-1 


1=^+1  *=-°°  1=^+1 
even  I  odd  I 

(4.151) 


4.5  Specific  Distributions 

Finding  actual  values  for  the  mean  and  variance  requires  a  specific  transformation  kernel. 
This  section  numerically  examines  the  mean  and  standard  deviation  obtained  when  using 
two  members  of  Cohen’s  class,  the  Wigner  distribution  and  the  binomial  distribution.  The 
analytic  results  derived  include  both  the  alias-free  and  the  traditional  formulations  of  these 
distributions. 


4,5- 1  Discrete  Wigner  kernel 

The  discrete  Wigner  kernel,  or  perhaps  more  properly  the  quasi- Wigner  kernel  [30],  has 
the  alias-free  form 


V^(n,0  =  Sn 

=  ^<5n  + ^<571-1 


(4.152) 

(4.153) 


which  makes  the  alias-free  variance 
var[Cj(n,p)] 


=  a 


2-00  2-00  .  v; 

E  E  ^-k^-k+  E  E  (2^-*^ 


1=^+1  k=-cx> 
even  I 


k=-oo 

odd  I 


(4.154) 


=  a' 


even  I 


1  + 


odd  I 


1 

2 


(4.155) 


(4.156) 
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As  mentioned  in  section  4.4.1,  the  traditional  formulation  sets  the  function  ip'{n,l)  to 
zero.  This  forces  all  the  rows  in  the  ambiguity  domain  associated  with  odd  lags  to  zero. 
When  applied  to  the  variance,  this  gives 


vai[Cg{n,p)]  =  S  S-kS-k  (4.157) 

1=^+1  k=-oo 
even  I 

=  53  1  (4.158) 

1=^+1 
even  I 

=  a"(|-l)  (4.159) 

4.5.2  Binomial  Kernel 


Kernels  which  generate  a  reduced  interference  distribution  (RID),  such  as  the  binomial 
and  the  spectrogram,  provide  lower  values  for  variance  in  the  time-frequency  domain  than 
the  Wigner.  The  RID  kernels  replace  the  uniform  ambiguity  domain  weighting  of  the 
Wigner  kernal  with  a  tapered  window  before  transforming  the  values  into  the  time-frequency 
domain.  This  provides  a  measure  of  smoothing  to  the  final  distribution  at  the  cost  of  some 
frequency  and  time  resolution. 

Determining  the  variance  for  the  binomial  kernel  is  more  involved  since  unlike  the 
Wigner,  the  kernel  does  not  have  a  simple  delta  function  form  in  the  autocorrelation  domain. 
Exploiting  the  symmetry  in  the  kernel  allows  expressing  the  variance  as 


vw:[Cg{n,p)]  = 


--1 


E  E  E  E  (4-i60) 

=  ^+1  *=-oo  l==f+i 

ven  I  odd  I 

Each  row  of  the  binomial  kernel  is  given  by  the  binomial  coefficients,  scaled  so  the  sum 
across  the  row  is  imity.  Both  ‘ip  and  xp'  have  the  same  form,  allowing  the  two  summations 


to  be  combined. 


r-l 


yaj:[Cg{n,p)]  =  a'^  Z)  E  [V’(fc,0]^ 


(4.161) 


i=^+lA:=-oo 


Row  I  of  the  kernel  has  1 4- 1  non-zero  elements.  These  elements  have  the  unsealed  values 
based  on  the  binomial  coefficients 


\m  )  -  m)\ 


(4.162) 
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with  m  =  0 . . .  Each  element  in  row  I  must  be  scaled  by  2  ^  to  force  the  sum  of  the  values 
across  the  row  to  equal  1.  Replacing  the  k  summation  with  a  summation  over  m  gives 


var[C,(n.p)]  = ^  (4.163) 

This  is  symmetric  about  I  =  0  reducing  the  variance  for  the  alias-free  form  to 

^  jp  ^ 

var[C,(n,i,)|=„‘  1+2^  ^  E  [„,(/:„),]'[  (4164) 

I— I  m—0  ^  ^ 

w  j 

In  the  traditional  form,  the  even  lag  rows  have  the  same  values,  but  the  odd  lag  rows 
contain  only  zeros.  Since  summing  across  the  zero  rows  is  not  necessary,  the  variance  can 
be  re-written  as 


v^[C,{n,p)]=a^h+2  ^  ^  ^ 


(4.165) 


Figure  4.1  shows  how  the  variance  of  the  Wigner  and  binomial  distributions  vary  as  the 


Wigner  Variance  Binomial  Variance 


Transform  Length  Transform  Length 

Figure  4.1:  Noise  Variance  vs  Transform  Length 


length  of  the  transform  increases. 

4.6  Numerical  Results 

The  numerical  values  were  derived  from  a  vector  containing  4096  independent,  normally 
distributed  point  with  mean  0  and  variance  1.  The  time-frequency  distribution  was  found  for 
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this  vector,  using  traditional  and  alias-free  formulations  with  varying  frequency  resolution. 

The  values  for  the  time- frequency  distribution  filled  a  matrix  (7^(n,p).  The  trial  mean 
was  given  by 

mean  ^  E[Cg]^  (4-166) 

n  p 

where  N  is  the  niunber  of  time  samples,  in  this  case  4096,  and  P  is  the  transform  length 
and  the  niunber  of  frequency  bins  which  ranged  from  4  to  256. 

The  variance  was  estimated  by  squaring  each  element  in  the  time-frequency  domain, 
finding  the  average  of  these  squared  values,  and  subtracting  the  square  of  the  mean  value. 
This  can  be  expressed  as 

var  =  ^  ^  ^  Cg{n,pf  -  E  [Cgf  (4.167) 

n  p 

Tables  4.1  and  4.2  show  how  well  the  analytically  obtained  values  for  the  distribution 


Window 

Length 

Mean 

ideal  trial 

Traditional 

Variance 

ideal  trial 

Alias- Free 

Variance 

ideal  trial 

4 

1 

0.9978 

1 

1.0461 

2 

1.9973 

8 

1 

0.9978 

3 

3.0689 

5 

7.0689 

16 

1 

0.9978 

7 

7.0730 

11 

11.0606 

32 

1 

0.9978 

15 

14.8377 

23 

22.7618 

64 

1 

0.9978 

31 

30.7465 

47 

46.5828 

128 

1 

0.9978 

63 

62.1522 

95 

93.6993 

256 

1 

0.9978 

127 

124.6105 

191 

187.3551 

Table  4.1:  Wigner  Statistics 


statistics  match  the  values  obtained  through  a  computer  trial.  These  tables  show  excellent 
agreement  between  the  theoretical  and  experimental  statistics,  with  only  two  trial  variance 
values  differing  from  the  expected  analytical  values  by  more  than  four  percent. 

These  results  show  a  growth  in  variance  similar  to  that  displayed  when  computing  a 
periodogram  [31].  A  more  complete  understanding  would  be  provided  by  considering  the 
case  of  a  desired  signal,  s{t),  added  to  the  noise  considered  in  this  chapter.  This  would 
allow  evaluation  of  the  performance  of  Cohen’s  class  transformation  as  a  signal  estimator 
in  terms  of  bias  and  asymptotic  performance. 
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Window 

Length 

Mean 

ideal  trial 

Traditional 

Variance 

ideal  trial 

Alias-Rree 

Variance 

ideal  trial 

4 

1 

0.9978 

1.0000 

1.0461 

2.0000 

1.9973 

8 

1 

0.9978 

1.7500 

1.8054 

3.3750 

3.4049 

16 

1 

0.9978 

2.7480 

2.8194 

5.2842 

5.3348 

32 

1 

0.9978 

4.1145 

4.1613 

7.9568 

7.9594 

64 

1 

0.9978 

6.0152 

6.0542 

11.7164 

11.6903 

128 

1 

0.9978 

8.6812 

8.7077 

17.0188 

16.9346 

256 

1 

0.9978 

12.4538 

12.3190 

24.5074 

24.1765 

Table  4.2:  Binomial  Statistics 


4.7  Summary 

Noise  driving  a  Cohen’s  class  transformation  provides  a  time-frequency  distribution 
with  characteristics  based  on  the  kernel  used  and  the  length  of  the  transform.  Comparing 
the  quasi-Wigner  and  the  binomial  kernels  demonstrated  a  significant  edge  in  performance 
for  the  binomial  as  the  frequency  resolution  (and  the  transform  length)  increased.  The 
quasi-Wigner  displayed  shghtly  worse  than  linear  growth  in  the  distribution  variance  as  the 
transform  length  increased,  0{N),  while  the  binomial’s  variance  was  slightly  greater  than 

o{Vn). 
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CHAPTER  5 


DYNAMIC  TARGETS 


Studies  of  long  wires  and  pairs  of  orbiting  spheres  demonstrate  the  usefulness  of  time- 
frequency  distributions  for  the  analysis  of  radar  signals.  They  served  as  simplified  models 
for  moving  blades  of  an  engine  propeller  or  a  helicopter  rotor.  In  general,  the  backscatter 
signal  was  difficult  to  interpret  in  either  the  time  or  the  frequency  domains.  Applying  the 
binomial  distribution,  a  discrete  time-frequency  distribution,  allows  clearly  associating  each 
sphere  with  its  corresponding  doppler  return.  The  binomial  distribution  provided  a  detailed 
view  of  the  target  dynamics,  opening  the  way  for  target  classification  and  identification.  The 
structures  and  details  available  in  the  time-frequency  domain  were  not  readily  exploited  in 
the  original  signal  representation. 

5.1  Introduction 

In  addition  to  looking  at  dispersive  targets,  time-frequency  techniques  are  also  directly 
applicable  to  dynamic  targets.  With  the  target  in  motion,  the  scattered  signal  will  have 
non-stationary  characteristics.  Focusing  on  these  changing  characteristics  allows  estimation 
of  parameters  relating  to  the  target’s  geometry  and  insight  into  the  target  dynamics. 

Radar  systems  combine  elements  of  electromagnetic  theory  with  signal  processing  to 
estimate  target  parameters.  A  rudimentary  system  may  only  try  to  determine  whether  or 
not  a  target  is  present  in  a  region.  More  elaborate  systems  provide  estimates  of  the  target’s 
location  and  velocity  and  direction.  Some  systems  extend  the  parameter  estimation  to 
provide  target  identification  [36]. 

Simple  radar  analysis  can  determine  position  and  radial  velocity.  The  target  location 
is  resolved  in  polar  coordinates  by  using  time  delay,  AT,  in  the  returned  signal  to  measme 
range  and  the  pointing  direction  of  a  narrow-beam  antenna  to  provide  the  necessary  angular 
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coordinates  (azimuth  and  elevation).  Transforming  the  returned  signal  into  the  frequency 
domain  gives  the  doppler  shift  imposed  on  the  signal  by  the  target.  Knowing  this  shift 
and  the  original  carrier  frequency  gives  the  radial  velocity  of  the  target  (how  quickly  it  is 
moving  directly  towards  or  away  from  the  radar). 

While  these  techniques  work  reasonably  well  for  a  simple  return,  multiple  targets,  or 
multiple  scatterers  on  a  single  target,  can  obscure  the  analysis.  Interaction  terms  make 
it  difficult  to  determine  the  contribution  from  individual  scatterers  or  scattering  modes. 
Advanced  methods  like  time-frequency  analysis  can  isolate  components  from  such  multiple 
returns. 


5.2  Background 


The  basic  motion  selected  for  this  work  was  rotation  about  a  point.  The  targets  con¬ 
sisted  of  wires  and  spheres.  The  objective  for  the  time-frequency  analysis  was  to  determine 
instantaneous  target  velocities  based  on  the  doppler  shift  of  the  return  signal.  For  a  carrier 
frequency  of  /c,  the  doppler  shift  for  a  point  target  moving  with  a  radial  velocity  of  v  is 


c 


(5.1) 


using  c  as  the  velocity  of  propagation. 

Figure  5.1  shows  the  basic  measurement  system  modeled  to  measure  these  doppler  shifts. 
The  target  was  illuminated  with  a  1  GHz  plane  wave,  then  the  scattered  signal  from  the 
target  was  converted  down  to  the  baseband  in-phase,  /,  and  quadrature,  Q,  channels. 

For  an  object  following  a  circular  path,  the  radial  velocity  will  follow  a  sinusoid  with 
the  maximum  directly  proportional  to  the  rotational  rate,  /rot,  and  the  distance  from  the 
center  of  rotation,  r.  The  maximum  radial  velocity  is  27Tfrot'^^  so  the  maximum  doppler 
shift  is 


fd„ 


^'^frotf  fc 


(5.2) 


The  rotational  rate  and  distance  were  selected  in  conjunction  with  the  sampling  frequency 
to  avoid  any  ambiguities  in  the  instantaneous  doppler  frequency. 


5.3  Wire  Targets 

Thin,  straight  pieces  of  wire  provided  a  relatively  simple  yet  potentially  important  target 
for  analysis.  The  insights  provided  from  a  single  wire  carry  over  into  multiple  wire  targets. 
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Figure  5.1:  Dynamic  Measurement  System  Model 


Rotating  wires  serve  as  a  simplified  representation  of  rotating  structmes  such  as  airplane 
propellers  and  helicopter  rotors. 

5.3.1  Single  Rotating  Wire 

The  first  set  of  dynamic  targets  examined  were  long,  perfectly  conducting  wires,  modeled 
using  the  Numerical  Electromagnetics  Code  (NEC)  [5].  These  simple  targets  illustrate  how 
time-frequency  distributions  can  be  applied  to  a  radar  scattering  problem. 

Since  the  wires  used  were  electrically  long,  the  backscatter  pattern  has  several  expected 
peaks.  In  addition  to  a  sharp  peak  when  the  wire  is  broadside  to  the  incident  field,  two 
other  broad  peaks  result  due  to  traveling  waves  set  up  on  the  wire.  These  peaks  appear 
at  [24] 


For  the  wires  used,  2.5  and  5  wavelengths  in  length,  these  scattering  wave  peaks  occur  at 
22.1  and  31.2  degrees. 

These  long  wires  demonstrated  a  second  scattering  property.  The  primary  scattering 
center  was  located  at  the  tip  of  the  wire  furthest  firom  the  radar.  As  the  wire  rotates,  the 
primary  doppler  contributor  will  follow  the  movement  of  the  outer  tip  of  the  wire. 
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Three  different  single  wire  cases  were  tested.  In  two  of  the  trials,  the  wire  was  rotated 
around  one  end.  In  the  third  trial,  the  rotational  center  was  placed  2.5  wavelengths  beyond 
the  end  of  the  wire.  This  allowed  the  end  of  the  2.5  wavelength  wire  to  reach  the  same 
radial  velocity,  and  therefore  generate  the  same  doppler  shift,  as  the  5  wavelength  wire. 
Figure  5.2  illustrates  these  cases. 


Figure  5.2:  Single  Wire  Geometries 


The  corresponding  backscatter  patterns  for  the  2.5  wavelength  and  the  5  wavelength 
wires  are  shown  in  figure  5.3 

2.5  Wavelsnglhs  5  Wavelengihs 
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Figure  5.4  shows  the  magnitude  of  the  time  response  for  these  three  cases.  A  significant 


2.5  Wavelength 


Figure  5.4:  Single  Wire  Waveforms 

feature  is  how  both  2.5  wavelength  targets  generate  the  same  received  signal  power  as  a 
function  of  time.  The  difference  in  rotational  centers  only  appears  in  the  phase  of  the 
received  signal. 

Taking  the  Fourier  transform  of  the  signals  provides  the  power  spectral  densities  shown 
in  figure  5.5.  Because  the  scattering  centers  on  the  wire  have  constantly  changing  velocities, 
the  doppler  of  the  scattered  signal  has  a  smeared  nature,  covering  a  range  of  values.  The 
differences  between  the  2.5  and  5  wavelength  wires  can  be  seen  in  the  maximum  doppler 
achieved,  with  the  5  wavelength  wire  having  twice  the  frequency  spread. 

Another  interesting  feature  in  figure  5.5  is  the  difference  between  the  2.5  wavelength 
wires.  As  expected,  the  wire  with  the  offset  rotational  center  reaches  the  same  maximum 
frequency  as  the  longer  wire.  Both  ends  of  the  shorter  wire  are  in  motion,  so  neither  tip 
stays  stationary  at  the  rotational  center  to  make  a  significant  contributions  at  the  lower 
doppler  frequencies. 

The  binomial  distribution  for  the  single  wires  is  shown  in  figure  5.6.  In  the  joint  time- 
frequency  domain  several  characteristics  of  the  signal  become  obvious.  The  rotational  aspect 
of  the  wire  dynamics  show  up  in  the  sinusoidal  envelope  formed  in  the  binomial  distribution. 
The  maximum  doppler  frequency  occurs  when  the  wire  is  perpendicular  to  the  direction  of 
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Figure  5.6:  Single  Wire  Binomial  Distributions 


'4 


propagation,  so  it  corresponds  to  one  of  the  peak  amplitudes  in  the  return.  The  broadside 
orientation  also  allows  each  point  along  the  wire  to  contribute  energy  to  the  return,  creating 
a  broadband  frequency  return  as  all  dopplers  between  0  and  fmax  contribute. 

The  largest  returns  in  the  backscatter  signal  are  due  to  the  traveling  wave  returns.  These 
show  up  four  times  in  each  cycle.  The  extreme  size  of  this  component  obscures  the  details 
of  the  time-frequency  distribution  around  these  points  in  addition  to  creating  cross  terms. 

Another  characteristic  of  the  scattered  wave  shown  by  the  time-frequency  distribution 
is  how  the  scattering  favors  the  end  of  the  wire  furthest  from  the  radar.  On  these  plots, 
the  sinusoidal  pattern  generated  is  asymmetric  through  a  cycle.  The  trace  is  darker  from 
the  negative  maximum  doppler  point  through  the  positive  maximum  doppler  one  half  cycle 
later.  This  corresponds  to  the  time  the  outer  tip  is  further  from  the  radar  than  the  inner 
point.  The  other  half  cycle  has  a  lighter  component  and  corresponds  to  the  time  when 
the  fast  moving  outer  point  is  closer  to  the  radar.  For  the  wires  rotated  about  their  end, 
the  zero  frequency  line  is  stronger  during  this  half  cycle.  For  the  offset  wire,  there  are 
no  significant  contributions  from  the  spectrum  around  zero  hertz.  However,  there  is  still  a 
jump  in  emphasis  following  the  rear  point.  In  this  case,  neither  end  of  the  wire  is  stationary, 
so  the  zero  doppler  contribution  is  very  small. 

The  offset  wire  has  an  interesting  overall  pattern.  It  has  the  same  basic  features  as  the 
other  2.5  wavelength  wire,  but  a  null  band  appears  at  the  low  doppler  frequencies.  At  the 
maximum  doppler  points,  the  offset  wire  reaches  the  same  peak  frequency  as  the  longer 
wire,  as  expected,  but  it  cannot  support  modes  with  zero  doppler  shifts  since  all  the  points 
on  the  offset  wire  are  traveling  fast  enough  to  generate  a  doppler  shift  of  half  the  peak  value. 

5.3.2  Multiple  Rotating  Wires 

The  scattering  off  a  2.5  wavelength  wire  was  extended  to  consider  multiple  wires.  In 
this  set  of  trials,  the  2.5  wavelength  section  served  as  the  initial  target.  Then  a  second, 
identical,  wire  was  attached  at  the  rotation  point,  positioned  120  degrees  away  from  the 
first  wire.  A  third  wire  attached  to  the  rotation  point  and  positioned  240  degrees  away 
from  the  first  wire  completed  the  symmetric  three-wire  set.  The  geometries  involved  are 
shown  in  figure  5.7.  The  backscatter  for  the  single  wire  was  already  shown  in  figure  5.3. 
The  corresponding  backscatter  patterns  for  the  double  and  triple  wire  targets  are  given  in 
figure  5.8. 

The  targets  were  rotated  while  illuminated  with  a  planar  incident  field  with  the  electric 
field  polarized  parallel  to  the  plane  of  the  wires.  Figure  5.9  shows  the  magnitude  of  the 
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Double  Wire  Target 
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Single  Wire 


Figure  5.9:  Multiple  Wire  Waveforms 


received  backscatter  signal.  The  corresponding  power  spectral  density,  obtained  using  the 
Fourier  transform,  is  shown  in  figure  5.10. 

As  with  the  single  wire  trials,  the  dominant  returns  are  due  to  traveling  waves.  With 
careful  processing  of  the  scattered  signal,  it  is  possible  to  discern  portions  of  two  sinusoidal 
patterns  in  the  two-wire  target.  This  is  clearest  in  the  vicinity  of  the  maximum  doppler 
shifts.  Figure  5.11  shows  how  increasing  the  complexity  of  the  target  return  rapidly  fills  the 
time-frequency  response  with  components  and  cross  terms  that  make  interpretation  difficult. 
The  rapid  changes  in  the  power  of  the  backscatter  signal  contribute  to  this  difficulty. 

There  is  one  additional  problem  encountered  when  trying  to  interpret  the  three-wire 
scatterer  compared  to  the  single  and  double  wire  cases.  The  three-wire  target  is  the  only 
one  in  the  set  that  exhibits  complete  rotational  symmetry  —  each  third  of  a  rotation  is 
indistinguishable  from  the  other  two  thirds.  This  makes  it  extremely  difficult  to  make 
analysis  without  a  priori  knowledge  of  the  target  makeup. 

5.4  Spherical  Targets 

This  study  used  two  spheres  traveling  in  simple  circular  orbits.  Although  the  spheres 
were  illuminated  by  a  planar,  continuous-wave  signal,  the  return  from  the  two  moving 
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Frequency  Index 


Figure  5.10:  Multiple  Wire  Spectral  Densities 


Single  Wire  TFD 


Figure  5.11:  Multiple  Wire  Binomial  Distributions 
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Figure  5.12:  Target  Dimensions 

scatterers  undergoes  rapid  changes  in  both  amplitude  and  phase. 

Using  the  scattered  signal  allows  determining  the  number  of  scatterers,  their  rotational 
frequency,  and  position  of  the  spheres  with  respect  to  the  center  of  rotation.  By  using 
time-frequency  techniques,  the  typical  estimates  of  these  values  had  errors  less  than  five 
percent. 

5.4.1  Electromagnetics  Problem 

The  basic  problem  examined  involved  two  metallic  spheres  as  sketched  in  figure  5.12. 
The  target  consisted  of  two  perfectly  conducting  spheres  with  their  centers  separated  by  one 
meter.  Each  sphere  had  a  radius  of  15  centimeters.  When  illuminated  by  a  1  GHz  incident 
plane  wave,  the  multiple  spheres  created  several  different  modes  of  reflected  electromagnetic 
energy.  These  modes  included  specular  reflections,  creeping  waves,  and  multiple  bounce 
reflections. 

The  original  data  for  this  experiment  was  generated  by  running  a  body  of  revolution 
simulation  using  two  spheres  as  targets.  The  returned  signal  values  were  based  on  monos¬ 
tatic,  far-field  calculations.  The  simulation  covered  a  frequency  range  of  1.0  to  7.55  GHz  and 
considered  observation  angles  over  a  span  of  90  degrees.  Viewed  in  polar  coordinates,  the 
spheres  were  placed  with  one  center  at  (0.0, 0.0, 0.5)  and  the  other  center  at  (0.0, 0.0,  —0.5). 
The  angle  6  ranged  from  0  to  90  degrees,  sampled  every  half  degree.  The  symmetry  of  the 
two  spheres  allowed  this  data  set  to  be  extended  to  cover  view  angles  firom  0  to  360  degrees 
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Figure  5.13:  Centers  of  Rotation 


with  respect  to  the  positive  z-axis,  measured  in  the  x-z  plane.  The  carrier  frequency  of 
1  GHz  insured  the  doppler  shifts  imposed  on  the  signal  by  the  target  rotation  would  be 
small  enough  to  avoid  undersampling  in  the  frequency  domain. 

5.4.2  Processing  and  Estimation 

To  create  a  dynamic  problem,  numerical  simulation  placed  the  spheres  in  circular  orbits 
around  five  different  rotational  centers,  referred  to  as  cases  A  through  E.  Figure  5.13  shows 
the  positions  chosen.  Each  of  the  five  cases  provided  different  return  signal  characteristics. 
Each  trial  used  two  seconds  of  data  sampled  at  a  rate  of  512  hertz.  After  applying  a  512- 
point  transformation,  this  provided  resolution  of  1  hertz  in  the  frequency  domain.  The 
rotational  rate  was  chosen  as  3.5  hertz  so  the  maximum  expected  doppler  shift  would  not 
give  an  undersampled  signal.  Calculations  for  each  sampling  instant  determined  the  view 
angle,  then  interpolated  between  the  values  in  the  simulated  data  set.  In  all  cases,  each 
sphere  maintained  a  constant  distance  from  the  rotational  center.  Both  spheres  traveled  at 
the  same  angular  rate. 

For  each  of  the  five  cases,  estimates  were  generated  of  four  parameters:  the  period 
(To),  the  distance  from  the  rotational  center  to  the  first  sphere  (ri),  the  distance  from  the 
rotational  center  to  the  second  sphere  (r2),  and  the  distance  between  the  spheres  (ro). 

The  period  was  determined  using  the  autocorrelation,  R(t),  of  the  time  domain  samples. 
The  number  of  autocorrelation  peaks,  N,  and  the  time  shifts  for  the  first  and  last  peaks,  ri 
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and  Tat,  provided  the  estimate  of  the  period 


_  Tf{  —  T\ 

N-\ 


(5.4) 


The  radii  of  the  sphere  orbits  were  based  on  the  measured  peak  doppler  shift,  fdmaxi 
the  estimated  rotational  period,  To,  and  the  known  illumination  frequency,  /«.  Using  equa¬ 
tion  5.2,  the  only  unknown  value  after  measuring  the  peak  doppler  shift  was  the  radius,  r. 
The  estimate  for  the  radius  rx  or  r2  becomes 


ri,2 


fdma.TpC 

47r/i 


(5.5) 


The  final  value,  ro,  was  found  using  ri,  r2,  and  the  phase  shift,  0,  between  their  peak 
doppler  points.  The  estimate  of  0  depends  on  Ti  and  T2,  the  times  each  sphere  generates 
their  maximum  doppler  shift 

0  =  SQollIlIl  (5.6) 

Using  this  phase  shift  in  the  law  of  cosines  gives  the  sphere  separation,  ro  as 


rp 


2rir2  COS0 


(5.7) 


When  generating  these  estimates,  the  period  was  calculated  autonomously  by  a  com¬ 
puter.  The  maximum  doppler  shift  frequencies  and  times  values  were  manually  extracted 
from  plots  of  the  time-frequency  distribution  to  allow  separating  the  doppler  components 
associated  with  each  sphere. 

Since  the  spheres  traveled  in  a  circular  path,  their  radial  velocity  with  respect  to  a 
stationary  observer  varied  sinusoidally.  The  velocity  passed  through  zero  at  the  spheres’ 
closest  and  furthest  position  with  respect  to  the  observer,  that  is,  when  the  sphere,  the 
rotational  center,  and  the  observer  were  in  line.  The  maximum  velocity  points  were  offset 
by  ±90  degrees  from  the  zero  velocity  points. 

The  doppler  shift  imposed  on  the  scattered  signal  is  directly  proportional  to  the  sphere 
velocity,  and  is  given  by  the  relationship 

/d  =  y  (5.8) 

Since  the  two  targets  cannot  occupy  the  same  position  of  the  same  orbit  simultaneously, 
each  will  have  a  different  doppler  frequency  over  the  course  of  a  cycle.  In  the  cases  we 
considered,  the  doppler  returns  will  all  have  the  same  period,  but  will  vary  in  their  peak 
firequency  deviations  and  their  relative  phases. 

Figure  5.14  shows  a  segment  of  the  backscatter  waveform.  Since  the  test  signals  are 
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Figure  5.14:  Orbiting  Spheres  Waveform 

distinguishable  only  by  their  phasing,  this  figure  represents  the  magnitude  of  the  complex 
value  signal  for  all  five  cases. 

The  test  signals  show  dramatic  differences  when  viewed  in  the  frequency  domain.  Two 
examples  of  the  fast  Fourier  transform  on  the  signal  gives  the  spectral  magnitude  shown  in 
figures  5.15  and  5.16.  Unfortunately,  the  continuously  changing  doppler  frequencies  cannot 
be  discerned  from  these  graphs.  These  spectral  graphs  can  also  used  to  estimate  the  period 
and  the  maximum  doppler  frequency  for  the  signal. 

Figure  5.17  shows  one  of  the  possible  spectrograms  for  the  Case  E  signal.  This  was 
calculated  using  a  128  point  window,  giving  a  frequency  resolution  of  4  hertz  (compared 
with  the  1  hertz  resolution  provided  by  a  512  point  binomial  transformation).  As  the  figure 
shows,  this  form  makes  it  difficult  to  associate  specific  frequencies  with  particular  times. 

The  binomial  distribution  gives  a  different  perspective  by  taking  the  signal  into  the 
time  frequency  domain,  shown  in  figures  5.18  and  5.19.  These  graphs  show  two  sinusoidal 
components,  one  created  by  each  sphere,  distinguished  by  their  magnitudes  and  phases. 
These  differences  translate  into  distance  firom  each  sphere  to  the  center,  and  from  the 
spheres  to  each  other. 
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Figure  5.15:  FFT  of  Case  A  Waveform 


Figure  5.16:  FFT  of  Case  E  Waveform 
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Figure  5.18:  Case  A  Time-Frequency  Representation 
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Figure  5.19:  Case  E  Time- Frequency  Representation 


To 

n 

r2 

^0 

Actual 

146.29 

0.50 

0.50 

1.00 

Estimate 

146.33 

0.50 

0.50 

1.02 

%  Error 

0.03 

0.98 

0.98 

1.97 

Table  5.1:  Case  A  Results 


5.4.3  Test  Results 
Case  A  —  Symmetric 

The  first  rotating  sphere  target  had  the  rotational  center  at  the  midpoint  between  the 
spheres,  indicated  by  position  A  in  figure  5.13.  This  required  no  additional  manipulation 
of  the  signal  since  the  rotational  center  matches  the  original  simulation. 

Since  the  spheres  were  equidistant  from  the  center  and  traveling  at  the  same  angular 
rate,  the  velocity  of  one  sphere  was  always  the  negative  of  the  other.  This  means  the  two 
spheres  impose  equal  and  opposite  doppler  shifts  on  the  returned  signal. 

Table  5.1  summarizes  the  estimation  results.  Due  to  the  high  level  of  symmetry  in  this 
case,  a  priori  knowledge  that  the  value  of  Tq  corresponded  to  half  the  actual  period  was 
used  when  generating  the  estimates. 
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To 

ri 

r2 

ro 

Actual 

146.29 

0.25 

0.75 

1.00 

Estimate 

146.25 

0.27 

0.76 

1.05 

%  Error 

-0.02 

6.44 

0.98 

4.75 

Table  5.2:  Case  B  Results 


To 

n 

r2 

ro 

Actual 

146.29 

0.00 

1.00 

1.00 

Estimate 

146.33 

0.00 

1.00 

1.00 

%  Error 

0.03 

0.00 

0.30 

0.30 

Table  5.3:  Case  C  Results 


5.4.4  On-axis  Shifts 

The  next  three  cases  share  a  common  difference  with  the  original  data  set.  Case  B  moved 
the  rotational  center  closer  to  one  sphere.  In  case  C  the  rotational  center  was  colocated 
with  the  center  of  one  sphere.  Case  D  positioned  the  rotational  center  on  the  extension  of 
the  line  connecting  the  sphere  centers.  In  these  three  cases,  the  return  signal  value  from 
each  angle  requires  an  additional  phase  factor,  to  account  for  the  apparent  movement  of 
the  original  rotational  center  and  the  wavelength.  A,  of  the  illuminating  signal.  To  move 
the  origin  p  meters  along  the  axis  between  the  spheres  used  a  phase  offset  a  given  by 

a  =  cos(0)  (5.9) 

The  larger  errors  in  table  5.2  are  due  to  difficulties  in  estimating  small  lengths  due  to 
the  small  doppler  shifts  they  produce.  This  same  effect  also  contributes  to  the  large  errors 
encountered  in  Case  E  below.  The  results  for  Case  C,  summarized  in  table  5.3,  and  Case  D, 
siunmarized  in  table  5.4,  had  lower  errors,  corresponding  to  the  larger  radii  being  estimated. 


To 

r\ 

r2 

ro 

Actual 

146.29 

0.50 

1.50 

1.00 

Estimate 

146.33 

0.50 

1.51 

1.02 

%  Error 

0.03 

-0.38 

0.53 

2.14 

Table  5.4:  Case  D  Results 
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To 

n 

r2 

ro 

Actual 

146.29 

0.38 

0.92 

1.00 

Estimate 

146.25 

0.41 

0.92 

1.13 

%  Error 

-0.02 

6.92 

-0.35 

12.58 

Table  5.5:  Case  E  Results 


5.4.5  OflF-axis  Shift 


The  rotational  center  is  placed  at  an  arbitrary  position,  specified  as  an  angle  and  distance 
from  the  symmetric  center.  This  case  requires  two  parameters  to  determine  the  phase 
corrections.  Like  the  previous  case,  the  distance  p  factors  in,  but  it  also  involves  the  angle 
Oof!  between  the  z-axis  and  the  line  connecting  the  original  and  new  centers  of  rotation. 
This  yields  a  phase  correction  of 


a  = 


^cos(0  +  0o//) 


(5.10) 


As  shown  in  table  5.5,  Case  E  exhibited  the  largest  errors,  in  part  because  it  was  the 
only  case  where  requiring  estimates  of  the  angular  offset  between  the  two  time-frequency 
domain  traces.  In  the  previous  four  cases,  we  place  the  doppler  components  either  exactly 
in  phase,  or  180  degrees  out  of  phase.  For  the  arbitrary  rotational  center  placement,  the 
error  estimating  the  phase  difference  combines  with  any  errors  estimating  ri  and  r2  to  give 
the  total  error  in  the  tq  value. 


5.4.6  Higher  Order  Effects 

There  are  several  higher  order  effects  that  can  arise  in  a  dual  sphere  problem.  The 
spherical  targets  will  support  a  creeping  wave,  where  part  of  the  impinging  energy  follows 
the  surface  of  the  sphere,  shedding  energy  as  it  travels  around,  some  of  it  directed  back 
towards  the  source.  Energy  following  this  path  can  be  distinguished  from  the  primary 
returns  due  to  an  additional  time  delay  from  their  longer  path.  They  do  not  show  up  with 
this  time-firequency  analysis  because  the  primary  (or  specular)  retmn  and  the  creeping  wave 
return  undergo  the  same  doppler  shift.  Note  that  at  extremely  high  sampling  rates,  two 
doppler  lines  could  arise,  offset  in  time  by  the  additional  travel  time  of  the  creeping  wave. 
For  a  15  centimeter  radius  sphere,  the  wave  travels  an  additional  0.15(7r-|-2)  =  0.77  meters, 
or  2.57  microseconds  in  free  space.  This  would  require  a  sampling  rate  over  388  MHz. 

Another  observable  effect  involved  waves  reflecting  off  both  spheres  in  succession  before 
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return  to  the  source.  Due  to  the  motion  of  the  reflection  points  on  the  surface  of  the  sphere, 
the  path  length  traveled  by  the  wave  between  the  spheres  varies  with  the  observation  angle, 
giving  a  roughly  hnear  value  over  each  half  cycle  of  a  rotation.  The  energy  involved  in  this 
mechanism  is  very  low  compared  to  the  specular  peaks,  and  does  not  show  up  in  most  of 
the  sample  runs. 

5.4.7  Customized  Kernel 

Since  this  problem  was  constrained  to  simple  rotational  motion,  it  is  possible  to  design 
a  kernel  optimized  to  extract  the  expected  doppler  shapes.  I  did  not  explore  that  option 
since  the  rotating  spheres  were  only  serving  as  an  example  to  justify  research  using  more 
complicated  targets.  Customizing  the  kernel  to  the  problem  may  provide  dramatic  improve¬ 
ments  in  the  resolution  and  fidelity  of  the  results,  particularly  when  the  signal  is  buried  in 
noise  [45]. 

5.5  Summary 

Time-frequency  techniques  worked  very  well  on  this  problem.  Examining  the  binomial 
distribution  allowed  us  to  extract  information  not  readily  available  in  other  representations 
of  the  signal.  FVom  this  information,  we  could  estimate  the  physical  parameters  of  the 
target  with  good  accuracy,  with  only  two  estimates  out  of  the  group  of  twenty  exceeding  a 
five  percent  error. 

The  techniques  used  here  are  applicable  to  other  types  of  problems.  Passive  systems 
which  measure  emitted  energy  could  also  use  this  type  of  time-frequency  analysis.  In  as¬ 
tronomy,  an  analogous  problem  might  be  the  signals  obtained  from  a  binary  star  system. 
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CHAPTER  6 


EFFECTS  OF  ADDING  NOISE  AND  COMPLEXITY  TO 

SCATTERED  SIGNALS 


The  previous  chapters  presented  the  application  of  time-frequency  techniques  to  several 
electromagnetic  problems.  While  the  results  have  been  promising,  their  applicability  is 
limited  due  to  the  idealized  conditions  employed.  Using  simulations  for  the  experiments 
provided  signals  free  of  the  noise  components  with  which  all  physical  systems  must  contend. 
Evaluating  this  aspect  is  essential  since  other  promising  techniques  for  scattered  signal 
analysis,  such  as  the  singularity  expansion  method  (SEM),  have  shown  a  loss  of  effectiveness 
once  a  realistic  level  of  noise  was  added  to  the  system  [4,  26]. 

For  this  chapter,  noise  signals  were  inserted  at  various  power  levels  to  determine  the 
robustness  of  the  time-frequency  techniques.  The  signal  degradation  was  monitored  using 
the  Renyi  third  order  entropy  [44]  and  the  effects  of  increasing  the  target’s  complexity 
examined. 


6.1  Background 

An  electronic  system,  especially  one  based  on  receiving  scattered  electromagnetic  waves, 
must  deal  with  a  variety  of  noise  sources.  These  sources  may  include  thermal  noise,  shot 
noise,  atmospheric  noise,  solar  noise,  cosmic  noise,  and  urban  noise.  In  some  instances, 
these  contributions  can  overpower  the  desired  signal,  making  it  impossible  to  extract  any 
useful  information.  While  some  systems  can  increase  their  signal  power  to  make  the  noise 
contribution  insignificant,  many  radars  cannot  do  this  due  to  size,  cost,  or  power  soiurce 
limitations.  Instead,  the  signal  processing  subsystems  in  the  radar  must  compensate  to 
extract  useful  information  from  the  imperfect  received  signal. 

To  investigate  the  impact  of  noise,  a  simplified  scattering  model  using  a  point  target 
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was  first  investigated.  This  scattered  signal  was  then  corrupted  by  a  noise  source,  in  this 
case  modeled  using  band-limited  additive  white  Gaussian  noise  (AWGN).  The  third-order 
Renyi  entropy  provided  a  quantitative  measure  of  the  noise’s  impact.  Various  amounts  of 
noise  were  then  added  to  corrupt  the  signal  to  see  how  the  received  information  degrades 
as  the  signal  to  noise  ratios  (SNR)  decreases.  Finally,  these  techniques  were  repeated  to 
corrupt  the  signals  from  the  two  orbiting  spheres  discussed  in  chapter  5. 


6.2  Signal  Models 

This  study  examined  four  target  sets,  two  based  on  ideal  scatterers  and  two  based 
on  realistic  physical  models.  The  data  sets  represent  the  scattered  field  from  a  target 
positioned  along  the  y-axis  and  rotated  around  the  z-axis.  The  targets  were  rotated  at  a 
rate  of  3  revolutions  per  second  and  illuminated  with  a  2  GHz  carrier.  The  scattered  field 
was  sampled  for  one  second  at  512  hertz.  Placing  the  target  pieces  within  0.75  meters  of 
the  origin  yields  a  maximum  expected  doi 

fd  = 


ipler  shift  of 


2v 

A 

(6.1) 

47r3(0.75) 

0.15 

(6.2) 

188.5  hertz 

(6.3) 

6.2.1  Point  Targets 


To  focus  on  the  noise  effects,  the  simplest  possible  data  signal  was  used.  Therefore, 
instead  of  the  accurately  modeled  sphere  scattering  discussed  earlier,  a  point  target  was 
placed  in  the  far  field  of  the  radar  antenna,  traveling  in  a  circular  orbit  of  radius  r  as 
illustrated  in  figure  6.1.  Since  the  target  is  in  the  far  field,  the  return  amplitude  remains 
constant  as  the  range  changes  slightly  with  target  motion.  The  target  motion  is  significant 
with  respect  to  illuminating  signal  wavelength  so  the  phase  will  depend  on  the  target-radar 
separation.  The  target’s  angular  position  at  time  t  depends  on  the  frequency  of  rotation 
ifrot)  and  is  given  by  2'irfrott- 

For  the  carrier  frequency  (/c),  the  signal  will  have  a  wavelength 

300,  OOOkm/s 


/c 

This  means  the  signal  will  undergo  a  round-trip  phase  shift  (0)  of 

^  47rcos(27r/roft) 


(6.4) 


(6.5) 
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Figure  6.1:  Point  Target  Geometry 


radians.  The  returned  signal  is 

git)  =  (6.6) 

With  the  target  traveling  in  a  circular  path,  the  range  to  the  target  is  given  by  a 
sinusoidal  function.  This  means  the  time-frequency  distribution,  based  on  the  velocity,  or 
the  derivative  of  the  range,  is  also  given  by  a  sinusoid.  The  doppler  shift  imposed  on  the 
signal  is  based  on  the  instantaneous  radial  velocity 

Vradit)  =  r  sm(27r/rott)  (6.7) 


and  the  carrier  frequency  wavelength 


fdop{i)  — 


Ac 


(6.8) 


Figure  6.2  shows  the  difference  between  a  traditional  Fourier  transform  and  a  time 
frequency  transformation  for  an  orbiting  point  target.  The  binomial  distribution  clearly 
shows  how  the  target’s  doppler  changes  as  a  function  of  time. 

The  first  two  targets  consisted  of  idealized  point  scatterers.  The  one  point  target  placed 
the  scatterer  at  the  point  y  =  0.75  (0.75  meters  from  the  center  of  rotation).  The  scattered 
signal  from  the  single  point  target  was  imique  since  it  exhibited  no  power  fluctuations,  only 
variations  in  phase.  The  two-point  target,  shown  in  figure  6.3,  used  this  scatterer  and  added 
a  second  scatterer  at  y  =  —0.25  (0.25  meters  on  the  opposite  side  of  the  center  of  rotation). 
Due  to  constructive  and  destructive  interference  between  the  scatter  fields  from  the  two 
points,  variations  occur  in  both  the  magnitude  and  phase  of  the  signals. 

The  simplified  nature  of  these  signals  made  them  suitable  for  evaluation  using  the  re¬ 
ceiver  operating  characteristic.  Calculating  the  receiver  operating  characteristic  requires 
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Frequency  (Hz)  Time  (sec) 

Figure  6.2:  Transformations  for  Single  Point  Scatterer 


.25  meter 

Figure  6.3:  Two-Point  Scattering  Geometry 
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.75  meter 


.15  meter 

Figure  6.4:  Two-Sphere  Scattering  Geometry 


an  apriori  time-frequency  distribution  to  determine  where  the  targets  components  should 
appear.  This  definitive  template  for  distinguishing  false  alarms  from  legitimate  detection  is 
only  possible  for  the  simplest  targets. 

6.2.2  Distributed  Target 

The  other  two  targets  were  more  difficult  to  analyze  since  the  ideal  point  scatterers 
were  replaced  by  numerical  simulations  modeling  physical  targets  which  generated  more 
complicated  scattered  fields.  The  first  realistic  target,  shown  in  figure  6.4,  consisted  of  two 
metallic  spheres,  each  0.15  meters  in  radius.  The  centers  of  the  spheres  were  positioned  at 
y  =  0.75  and  y  —  —0.25  meters.  The  final  target,  shown  in  figure  6.5,  was  a  1  meter  long  wire 
placed  along  the  y-axis  between  the  points  y  =  0.75  and  y  =  —0.25  meters.  The  scattered 
signals  from  these  targets  were  not  suitable  for  evaluation  using  the  receiver  operating 
characteristic  since  one  could  not  dictate  a  priori  what  the  time  frequency  representations 
should  be  for  these  returns. 

Each  target  signal  was  normalized  to  have  a  mean  square  value  of  1.0.  To  create  a 
noise-corrupted  version,  the  target  signal  was  added  to  an  appropriately  scaled  noise  signal. 
The  noise  signal  was  modeled  as  complex  valued,  with  both  the  in-phase  and  quadrature 
channels  having  imcorrelated,  normally  distributed  values  with  mean  values  of  zero  and 
equal  variances. 
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.75  meter 


Figure  6.5:  Wire  Scattering  Geometry 


6.3  Noise  Model 

Adding  noise  to  the  signal  required  a  computational  method  of  modeling  the  noise 
sources.  My  goal  was  to  approximate  the  type  of  signal  degradation  that  would  be  found 
in  a  physical  radar  system.  The  first  step  in  developing  the  noise  model  was  to  identify 
the  expected  noise  sources  and  their  impact  on  the  received  signal.  The  second  step  was  to 
determine  a  mathematical  model  for  simulating  these  sources  of  noise  to  add  to  the  desired 
signals. 

6.3-1  Noise  Sources 

All  electrical  systems  are  subject  to  thermal  noise  due  to  the  random  movements  of 
electrons  within  the  system  components.  These  random  movements  are  directly  related 
to  the  temperature  of  the  components  since  temperature  on  the  macroscopic  level  is  due 
to  kinetic  energy  on  the  atomic  level.  Any  lossy  component  (including  the  antenna)  will 
generate  small  noise  signals  at  temperatures  above  absolute  zero.  A  special  case  of  thermal 
noise  called  shot  noise  affects  systems  which  deal  with  small  signals  and  large  amounts  of 
amplification  (like  a  radar).  This  comes  from  the  random  creation  and  loss  of  carriers  in 
the  amplifying  components  of  the  radar.  In  addition  to  its  own  thermal  noise,  an  antenna 
can  also  receive  external  noise  signals,  such  as  atmospheric  noise,  solar  noise,  cosmic  noise 
and  urban  (man-made)  noise.  In  most  instances,  thermal  noise  has  the  greatest  impact  on 
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system  performance  [27]. 


6.3.2  Additive  White  Gaussian  Noise 

While  removing  noise  from  the  random  movements  of  electrons  and  charge  carriers  is 
impossible,  modeling  their  impact  on  a  system  employs  a  straight-forward  technique.  Con¬ 
cern  here  was  with  the  combined  effect  of  an  uncounted  number  of  individual  contributors, 
all  following  the  same  statistical  model.  Combining  all  these  independent  contributors  and 
applying  the  central  limit  theorem  leads  to  an  approximately  gaussian  distribution  for  the 
noise  voltage. 

Since  the  signals  of  interest  are  complex  valued,  the  noise  model  must  provide  complex 
valued  samples.  The  received  signal  can  be  viewed  as  two  orthogonal  channels:  the  in-phase 
channel,  /,  corresponding  to  the  real  portion  of  the  complex  signal;  and  the  quadrature 
channel,  Q,  corresponding  to  the  imaginary  portion  of  the  complex  signal.  For  a  narrow- 
band  additive  white  Gaussian  noise  source,  the  noise  signal  is  given  by 

n{t)  =  nc{t)  cos  27r/ot  -I-  ns{t)  sin  27r/ot  (6.9) 

where  /o  is  the  operating  frequency  of  the  radar,  and  the  functions  ndt)  and  ns{t)  repre¬ 
sent  a  Gaussian  distributed,  zero-mean  random  amplitude  for  the  in-phase  and  quadrature 
channels.  With  orthogonal  I  and  Q  channels,  the  values  ndt)  and  ns{t)  are  independent 
and  identically  distributed.  The  noise  power  is  given  by 

Pn  =  E[n\t)]  =  E[nl{t)]  +  E[n1{t)]  =  2a^  (6.10) 

While  the  noise  in  each  channel  has  a  Gaussian  distribution,  the  magnitude  of  the 
noise,  \n{t)\  =  \/|n^P"+"in7p,  follows  a  Rayleigh  distribution,  and  the  phase  is  uniformly 
distributed  between  0  and  2it  radians. 

This  representation  of  additive  white  noise  has  assumed  a  band-limited  receiver.  The 
receiver  band-limits  the  noise  using  a  low-pass  filter  after  mixing  the  scattered  return  with 
the  carrier  frequency. 


6.4  Renyi  Entropy 

With  all  the  processing  performed  on  the  scattered  signals,  one  needed  a  method  to  eval¬ 
uate  the  resulting  time-frequency  distributions.  At  first  viewed  the  processed  time-frequency 
distributions  were  studied  and  qualitatively  judged  how  well  expected  features  appeared  in 
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the  processed  return.  While  this  qualitative  approach  could  distinguish  broad  differences 
when  processing  simple  signals,  it  was  inadequate  for  differentiating  the  noise  corrupted 
cases.  These  cases  required  a  more  objective,  quantitative  measure  of  performance.  This 
figure  of  merit  would  also  allow  the  development  of  automatic,  iterative  techniques. 

One  measure  of  the  complexity  of  a  time-frequency  distribution  is  the  entropy.  For  a 
two-dimensional  signal  Cs{t,f),  the  Shannon  entropy  [17]  would  be 

H{Cs)  =  -jJ  Cs{t,f)log2Cs{tJ)dtdf  (6.11) 


There  are  some  major  problems  with  this  entropy  measure.  Due  to  the  log2  operation, 
this  definition  will  only  apply  if  the  distribution  Cs  is  greater  than  zero  for  all  values  of  t 
and  /.  Zero  values  often  occur  in  time-frequency  distributions  due  to  the  limited  support 
of  time-frequency  components,  and  negative  values  are  a  result  in  many  distributions  which 
allow  negative  values  in  order  to  satisfy  the  marginal  distributions. 

To  avoid  zero  and  negative  values  in  the  entropy  calculation,  a  more  generalized  entropy 
definition  was  used.  For  a  bivariate  density,  P{x,y),  the  Renyi  entropy  is  given  by 


Ha{P)  = 


1  ,  S  S  P°'ix,y)dxdy 

-«  SSP{x,y)dxdy 


(6.12) 


When  a  =  1,  the  Renyi  entropy  matches  the  Shannon  definition.  By  allowing  the 
summation  of  values  to  occur  before  taking  the  logarithm,  the  Renyi  entropy  can  avoid  the 
zero  and  negative  values  that  cause  the  logarithm  to  fail. 

Converting  the  double  integrals  with  respect  to  time  and  frequency  (f  f  dxdy)  into  a 
double  summation  yields  the  Renyi  entropy  for  a  discrete  distribution 

fe])  =  ^  log, E  E  ( t'l)) 


This  form  is  easily  coded  in  a  computer.  To  compare  the  relative  entropy  between  two 
distributions,  the  log2  StSf  term  is  constant  and  can  be  dropped. 

The  Renyi  entropy  exhibits  five  key  properties  when  analyzing  multi-component  sig¬ 
nals  [3]: 


•  Ha  attempts  to  measure  the  number  of  components 


•  Ha  does  not  count  cross  components  for  odd  values  of  a  >  1 


•  Ha  has  an  upper  and  lower  bound 

•  Ha  values  are  invariant  to  signal  time  and  frequency  shifts 


96 


•  Ha  exhibits  extreme  sensitivity  to  the  phases  of  closely  spaced  components. 

The  first  four  properties  are  beneficial  for  this  application,  but  the  last  property  can 
cause  some  problems.  This  phase  sensitivity  manifests  itself  as  misleadingly  high  or  low 
entropy  values  for  a  distribution  due  to  the  contribution  of  cross-term  components.  Using 
reduced  interference  distributions  reduce  the  cross  terms  and  help  avoid  this  sensitivity. 

In  general,  a  can  take  on  any  positive  value  greater  than  or  equal  to  one.  This  restriction 
guarantees  the  existence  of  the  entropy  measure  [3].  However,  for  non-integer  values  of  a, 
the  entropy  will  yield  complex  numbers  of  limited  use.  When  a  takes  on  an  even  integer 
value,  the  entropy  expression  loses  the  ability  to  distinguish  positive  and  negative  values 
in  the  time-frequency  distribution.  This  has  led  to  the  suggestion  of  using  a  =  3  for  the 
entropy  calculations  [44]. 

In  this  study,  the  Renyi  entropy  values  provided  two  different  types  of  information. 
First,  the  entropy  represents  how  much  uncertainty  exists  in  the  distribution.  Noise  will 
increase  the  uncertainty,  so  a  good  processing  technique  will  reduce  the  noise  components 
in  the  final  distribution,  and  therefore  reduce  the  entropy.  The  second  use  for  the  Renyi 
value  is  to  identify  the  number  of  significant  signal  components.  The  larger  the  entropy 
value,  the  more  scatterers  can  be  placed  in  the  target  field.  With  the  proper  selection  of  6t 
and  Sf,  the  number  of  scatterers  can  be  estimated  using  [3] 

Nscat  =  (6.14) 

Unfortunately  this  value  is  valid  only  for  idealized  time-frequency  distributions.  Typical 
results  have  the  components  of  the  target  signal  spread  in  both  the  time  and  frequency 
directions,  resulting  in  a  reduction  of  the  signal  uncertainty,  and  an  underestimated  value 
for  iVgcat* 

6.4.1  Distribution  Evaluations 

After  finding  the  time  frequency  distribution,  a  numerical  metric  objectively  evaluated 
the  result.  The  first  measure  used  was  the  receiver  operating  characteristic  [38].  This 
involved  generating  an  ideal  distribution  for  comparison  with  the  signal’s  time-frequency 
distribution.  Applying  a  threshold  to  the  realized  signal  should  give  a  matching  distribution. 
Errors  can  be  expressed  in  radar  terms,  with  a  probability  of  detection  (Pa)  and  probability 
of  false  alarm  {Pfa)  for  each  threshold.  Plotting  P/a  versus  Pd  provides  a  way  to  compare 
different  distributions  and  the  impact  of  adding  different  amounts  of  noise.  Ideally,  the 
curve  would  allow  a  Pd  approaching  1.0  while  keeping  P/o  approximately  0.  Graphically, 
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the  closer  the  curve  gets  to  the  upper  left  corner  of  the  plot,  the  better  the  performance. 
Figure  6.6  shows  some  representative  receiver  operating  characteristic  curves. 

The  other  measure  used  was  the  Renyi  entropy  of  the  distribution  [3].  This  metric 
attempts  to  measure  a  distribution’s  information  content.  For  a  discrete  time  frequency 
distribution,  Cg[n,Tn],  the  Renyi  entropy  is  given  by 


m])  =  -i  log.  E  E  f  V 


En'  Em'Cg[n>,m'] 


3 

+  l0g2^t^/ 


(6.15) 


where  St  and  are  the  time  and  frequency  resolution. 

Since  entropy  measures  the  amount  of  disorder  in  a  system,  the  maximum  entropy  oc¬ 
curs  for  the  distribution  of  a  noise  only  signal.  As  the  signal  to  noise  ratio  improves,  the 
entropy  decreases.  Using  different  time-frequency  transform  kernels  generated  distribution 
with  widely  varying  entropy  values.  This  made  it  difficult  to  compare  values  from  different 
distributions.  For  example,  the  high  resolution  nature  of  the  Wigner  distribution  consis¬ 
tently  provided  lower  entropy  values  than  the  binomial  distribution  for  the  same  signal. 
Instead  of  the  absolute  entropy  value  from  different  transforms,  the  entropy  margin,  or  the 
entropy  improvement  (in  bits)  over  a  pure  noise  signal  was  emphasised. 


6.5  Analysis  and  Results 

The  analysis  for  each  target  considers  the  scatter  signals  with  signal  to  noise  ratios 
(SNRs)  between  10  down  to  —10  decibels,  with  the  special  cases  of  zero  noise  (SNR  =  oo) 
and  zero  signal  (SNR  =  —  oo).  They  were  transformed  to  the  time-frequency  domain  via 
the  short-time  Fourier  transform  (spectrogram),  the  binomial  transform,  and  the  Wigner 
transform. 


6.5.1  Receiver  Operating  Characteristics 

Figure  6.6  shows  some  representative  receiver  operating  characteristic  (ROC)  cmves. 
Each  plot  shows  the  ROC  for  the  spectrogram  (lowest  curve  in  each  set),  the  Wigner 
distribution  (center  curve  with  circled  points)  and  the  binomial  distribution  (top  curve). 
The  binomial  and  Wigner  distributions  significantly  outperformed  the  spectrogram,  with 
the  binomial  typically  holding  a  slight  edge.  The  differences  were  most  apparent  with  the 
two-point  scatterer  since  the  additional  signal  component  led  to  more  cross  terms  (and  more 
errors)  in  the  Wigner  distribution. 
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Two  Point  Tcurgot 


Figure  6.6:  Receiver  Operating  Characteristics 


One  Point 


Wire 


Figure  6.7:  Entropy  Margins 


6.5.2  Entropy  Measures 

The  performance  advantage  of  the  binomial  distribution  over  the  Wigner  for  complicated 
signals  was  also  apparent  from  the  entropy  calculations.  Figure  6.7  shows  the  entropy 
improvement  over  pure  noise  for  the  simple  single  point  scatterer  and  for  the  wire. 

With  the  single  point,  the  Wigner  (shown  with  circled  points)  was  the  top  performer, 
beating  the  binomial  at  high  SNRs.  But  as  the  SNR  decreased,  the  Wigner’s  performance 
dropped  to  roughly  match  the  binomial.  As  the  SNR  dropped  below  —3  dB,  the  smoothing 
properties  of  the  spectrogram  gave  it  a  slight  advantage,  although  all  three  distributions 
performed  poorly. 

For  the  complicated  scattered  signal  off  the  wire,  the  results  were  a  little  different. 
Once  again  the  spectrogram  yielded  the  smallest  entropy  improvement  over  pure  noise. 
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The  existence  of  significant  interfering  cross  terms  dropped  the  Wigner’s  entropy  margin 
(marked  with  circles)  almost  down  to  the  spectrogram  level.  The  binomial,  a  reduced 
interference  distribution,  performed  significantly  better  than  either  the  spectrogram  or  the 
Wigner.  This  advantage  was  evident  until  the  SNR  dropped  to  —10  dB,  at  which  point  the 
target  return  was  indistinguishable  from  the  noise. 

6.6  Summary 

The  characteristics  of  the  signal,  as  well  as  the  signal’s  quality,  as  measured  by  the  signal 
to  noise  ratio,  have  a  significant  impact  on  how  well  different  time-frequency  distributions 
will  perform.  The  receiver  operating  characteristic  can  give  a  clear  idea  of  the  transfor¬ 
mation  performance  but  requires  an  idealized  distribution  for  comparison.  The  inability 
to  determine  an  ideal  distribution  makes  the  receiver  operating  characteristic  difficult  to 
apply  except  for  simple  canonical  targets.  The  Renyi  entropy  values  did  not  require  any 
additional  information  about  the  signal  and  could  be  applied  to  any  signal,  giving  results 
consistent  with  qualitative  evaluations  of  the  distributions. 

The  work  presented  in  this  chapter  also  demonstrated  how  a  reduced  interference  dis¬ 
tribution  like  the  binomial  can  outperform  the  Wigner.  Although  the  Wigner  can  provide 
higher  resolution  with  respect  to  both  time  and  frequency,  this  resolution  comes  at  the 
expense  of  large  interfering  cross  terms  that  can  degrade  the  time  frequency  representation. 
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CHAPTER  7 


DEVELOPMENT  AND  USE  OF  A  CUSTOMIZED 

KERNEL 


The  Wigner  and  binomial  distributions  represent  general  purpose  time-frequency  ker¬ 
nels.  Although  they  may  provide  excellent  time-frequency  representations  for  many  signals, 
they  make  no  effort  to  exploit  characteristics  such  as  periodicity  which  may  exist  in  a  signal. 
This  section  examines  the  use  of  customized  time-frequency  kernels  which  improve  on  the 
performance  of  the  general  kernels. 

To  examine  how  customized  kernels  might  be  applied,  a  genetic  algorithm  was  developed 
to  maximize  the  Renyi  entropy  of  the  time-frequency  distribution.  The  kernel  design  was 
constrained  to  provide  a  kernel  with  reasonable  properties,  as  outlined  in  chapter  2. 

7.1  Background 

Due  to  the  generality  of  Cohen’s  class  of  transformation,  essentially  an  infinite  number  of 
different  time  frequency  representations  can  exist  for  a  given  signal  [11].  Based  on  the  results 
in  chapter  6,  the  goal  was  to  find  a  kernel  which  provided  time-frequency  distributions  with 
higher  Renyi  entropy  values  than  the  Wigner  or  the  binomial. 

7.1.1  Genetic  Algorithms 

Genetic  algorithms  are  an  optimization  technique  that  model  the  process  used  in  evolu¬ 
tion  and  natural  selection  [20].  Parameters  related  to  the  value  to  be  optimized  are  coded 
into  genes,  and  sets  of  genes  are  put  together  to  make  chromosomes.  Allowing  the  computer 
to  use  survival  of  the  fittest  to  update  the  available  gene  pool  serves  as  an  optimization 
technique.  This  has  been  applied  to  a  wide  variety  of  problems,  ranging  from  economic 
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prediction  [20]  to  pattern  classification  [35]  to  array  antenna  design  [16] » 

The  genetic  algorithm  was  chosen  for  this  problem  due  to  the  complicated  mapping 
between  the  parameters  driving  the  kernel  and  the  resulting  cost  function.  It  also  provided 
an  opportunity  to  apply  an  optimization  technique  that  has  not  been  used  in  this  area. 

Basic  Procedure 

After  deciding  how  to  encode  the  necessary  parameters  genetically,  the  procedure  fol¬ 
lowed  is  simple.  An  initial  random  population  is  generated.  These  genes  are  rank  ordered 
based  on  their  fitness  to  solve  the  selected  problem.  Those  scoring  the  highest  are  carried 
into  the  next  generation  or  iteration.  The  poorest  performers  are  discarded. 

The  next  step  is  to  fill  the  vacancies  in  the  population.  This  is  done  by  performing 
a  genetic  crossover  operation  between  two  parents  selected  from  the  previous  generation’s 
fittest  chromosomes. 

Once  this  new  generation  is  filled,  the  fitness  of  the  new  members  is  evaluated,  the 
chromosomes  are  rank  ordered,  and  the  process  of  discarding  members  and  creating  a  new 
generation  are  repeated. 

Terminating  the  Algorithm 

The  genetic  algorithm  is  a  global  search  technique  that  typically  performs  well,  but 
it  does  not  have  guaranteed  convergence.  The  code  implemented  used  two  criteria  for 
stopping.  The  first  was  based  on  the  top  performer  not  changing  over  a  given  number  of 
generations.  The  second  condition  limited  the  total  number  of  iterations. 

7.1.2  Problem  Design 

The  approach  followed  to  select  the  elementary  function,  h[t)^  needed  for  the  kernel  was 
to  use  a  genetic  algorithm  to  select  candidates.  According  to  the  design  rules  [22],  h{t) 
should  be  a  smooth  function  defined  over  the  range  —0.5  <  t  <  0.5  which  goes  smoothly  to 
zero  at  the  ends  when  |i|  =  0.5  and  has  an  area  of  1. 

The  elementary  function  was  defined  using  a  cosine  based  series 

N 

Ki)  =  H  «n  COS  ((2n  —  l)27ri)  (7.1) 

n=l 

where  there  are  N  weight  values,  Unj  which  are  assigned  values  in  the  range  0  <  o„  <  1. 
These  weights  axe  applied  to  the  odd  harmonics  of  the  cosine,  so  each  term  will  have  a 
maximum  at  t  =  0  and  go  to  zero  at  |t|  =  0.5, 
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The  target  used  for  these  trials  was  a  point  target  traveling  in  a  circular  path,  similar 
to  the  case  shown  in  figure  6.1.  The  point  was  rotated  at  5  hertz  a  distance  of  1  meter  from 
the  center  of  rotation  while  illuminated  with  a  1  gigahertz  plane  wave.  The  scattered  signal 
was  sampled  at  a  rate  of  512  hertz.  These  values  were  selected  to  provide  several  cycles 
over  a  1  second  record  without  doppler  frequency  ambiguities. 

7.2  System  Optimization  Results 

Once  the  basic  encoding  was  selected,  trials  were  conducted  to  see  how  well  the  genetic 
algorithm  could  select  a  kerneL  The  parameters  chosen  were  5  weighting  coefficients,  each 
coded  using  5  bits.  The  population  size  was  set  to  30,  with  a  maximum  of  15  generations 
before  terminating  the  genetic  algorithm.  While  using  more  weights,  more  bits,  a  larger 
population,  and  a  longer  running  algorithm  should  allow  greater  performance,  the  scope 
was  selected  to  show  the  feasibility  and  to  keep  the  problem  scaled  for  execution  on  a  486 
computer  system. 

To  generate  a  time- frequency  distribution  for  the  test  signal,  the  genetic  algorithm 
selected  the  elementary  function  shown  in  figure  7.1.  This  function  was  selected  over  a 


Elementary  Function 


Figure  7.1:  Genetically  Derived  Elementary  Function 


period  of  15  generations,  during  which  time  the  third-order  Renyi  entropy  of  the  resulting 
time-frequency  distribution  increased  from  an  initial  value  of  11.2000  to  the  final  value  of 
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11.3174.  For  comparison,  the  third  order  Renyi  values  for  the  Wigner  and  the  binomial 
distributions  of  this  signal  are  10.1904  and  11.0409. 

Figure  7.2  shows  the  resulting  time-frequency  distribution  for  the  point  scatterer.  This 


Genetic  Derived  1TD 


Figure  7.2:  Genetically  Derived  Time- Frequency  Distribution 

is  very  similar  to  the  distribution  obtained  using  the  binomial  kernel.  While  the  genetically 
derived  distribution  shares  the  desirable  distribution  properties  of  the  binomial  discussed  in 
chapter  3,  the  genetic  algorithm  did  provide  slightly  lower  cross  terms  near  the  maximum 
doppler  points. 

7.3  Summary 

The  results  of  using  the  genetic  algorithm  for  this  problem  were  very  promising.  The 
results  were  very  close  to  the  binomial  distribution,  showing  the  technique  could  meet  the 
standard  results.  This  approach  needs  to  be  taken  further  to  try  and  exploit  characteristics 
such  as  periodicities  in  the  signal.  It  should  be  possible  to  have  the  kernel  tuned  to  certain 
signal  parameters  such  as  the  general  shape  and  period  of  the  signal.  Used  with  a  genetic  al¬ 
gorithm,  this  would  provide  an  adaptive  processing  technique  giving  clearer  time-frequency 
distributions  than  the  standard  Wigner  and  binomial  kernels. 
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CHAPTER  8 


CONCLUSIONS 


This  work  was  interdisciplinary,  mixing  elements  and  concepts  from  electromagnetics 
with  those  of  signal  processing.  The  purpose  was  to  find  some  new  applications  for  some 
popular  and  powerful  signal  processing  techniques. 

8.1  Stationary  Targets 

Examining  stationary  targets,  using  a  sweep  of  frequencies,  allowed  estimating  the  range 
to  scattering  centers  on  the  target.  By  restating  the  transformations  used  in  time-fi:equency 
analysis  for  use  in  frequency-time  analysis,  it  was  easier  to  determine  what  scattering  mech¬ 
anisms  were  contributing  to  the  scattered  signal.  The  transformations  not  only  returned 
information  on  the  position  of  the  scattering  centers,  but  also  provided  insight  on  the  dis¬ 
persion,  if  any,  exhibited  by  the  scattering  mode. 

The  imderstanding  provided  through  this  type  of  analysis  can  be  applied  to  help  control 
the  scattering  from  a  target  in  either  a  monostatic  or  bistatic  situation. 

8.2  Dynamic  Targets 

The  situation  with  analyzing  a  dynamic  target  was  very  different.  Here  the  radar  gener¬ 
ated  a  continuous  wave  signal.  The  information  was  contained  in  the  doppler  shift  imposed 
on  the  signal  when  scattered  by  the  target.  Performing  time-frequency  analysis  on  the  re¬ 
ceived  signal  generated  a  distribution  showing  the  target’s  doppler  components  as  a  function 
of  time.  This  doppler  history  can  be  used  to  estimate  different  target  parameters. 
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8.3  Noise  Statistics 


Looking  at  how  noise  is  processed  by  a  transformation  belonging  to  Cohen’s  class  pro¬ 
vided  an  imderstanding  on  the  impact  of  noise  and  how  the  noise  in  the  time-frequency 
distribution  relates  to  the  transform  selected.  The  results  show  a  significant  advantage  for 
the  reduced  interference  distributions,  such  as  the  binomial,  over  the  standard  Wigner. 


8.4  Effects  of  Noise  and  Complexity 

Time-frequency  distributions  are  easiest  to  interpret  when  there  are  few  signal  com¬ 
ponents  and  a  very  high  signal  to  noise  ratio.  As  more  components  appear  in  the  time- 
frequency  distribution,  they  allow  the  formation  of  cross  terms,  typically  degrading  the 
clarity  of  the  distribution.  This  degradation  occurs  whether  the  additional  component  was 
actually  generated  by  the  target  or  was  simply  caused  by  noise. 

The  third  order  Renyi  entropy  proved  to  be  a  valuable  metric  for  numerically  evaluating 
the  quality  of  a  time-frequency  distribution.  These  entropy  values  corresponded  to  the 
subjective  quality  ratings  an  experienced  viewer  would  assign  to  the  different  distributions. 

The  results  showed  the  superiority  of  the  reduced  interference  distribution  except  when 
dealing  with  very  simple  targets  with  almost  no  noise  corrupting  the  signal. 

8.5  Customized  Kernels 

The  work  done  using  a  genetic  algorithm  provided  performance  exceeding  the  Wigner 
and  binomial  distribution  for  a  simple  target.  Before  exploring  this  area  more  thoroughly, 
careful  consideration  must  be  given  to  the  fitness  function  and  the  genetic  encoding  used. 

8.6  Value  of  Time-Frequency  Techniques 

This  research  has  shown  time-frequency  analysis  as  a  valuable  tool  for  analysis  of  radar 
and  other  electromagnetic  based  signals.  The  ability  to  identify  and  evaluate  different 
scattering  modes  provides  insight  not  readily  apparent  in  other  signal  representation. 

8.6.1  Image  Processing 

Once  a  signal  has  been  processed  using  a  time-frequency  transformation,  the  resulting 
function  of  two  variables  is  similar  in  many  aspects  to  an  image.  It  appears  that  applying 
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image  processing  techniques,  such  as  noise  reduction  or  edge  enhancement,  may  provide  the 
next  logical  step  for  processing  scattered  signals.  The  image  processing  approach  also  leads 
into  the  area  of  target  classification. 

8.6.2  Ultrawideband  Radar 

The  processing  of  ultrawideband  radar  signals  is  a  natural  extension  to  the  work  in 
this  dissertation.  Since  an  ultrawideband  radar  cannot  apply  narrowband  assumptions,  the 
ability  to  process  a  signal  in  a  joint  time-frequency  domain  should  prove  beneficial. 
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APPENDIX  A 


Notational  Conventions 


Continuous 

Discrete 

0(i/,  r) 

4>{q,l) 

ambiguity  domain  kernel 

ipit,  t) 

autocorrelation  domain  kernel 

Cg{t,f) 

Cg{n,p) 

time  frequency  distribution 

9{t) 

9{n) 

time  domain  signal 

9r(t) 

gr{n) 

real  portion  of  signal 

9iit) 

9i{n) 

imaginary  portion  of  signal 

G{f) 

Giq) 

frequency  domain  signal 

f 

9 

frequency  of  interest 

u,  A 

q,b 

dummy  frequency  variable 

t 

n 

time  of  interest 

u,v 

m,  c 

dummy  time  variables 

T,S 

Ud 

lag  variables 

<8>t 

E[x] 

^tf 


convolution  along  the  t  axis 
expected  value  of  x  (also  written  x) 

Fourier  transform  from  g{t)  to  G{f) 
inverse  Fourier  transform  from  G{f)  to  g{t) 
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APPENDIX  B 


Computer  Subroutines 


This  appendix  lists  the  major  Matlab  routines  used  in  this  research.  The  routines  assume 
the  signal  will  be  provided  as  a  one-dimensional  vector.  The  values  may  be  real  or  complex. 
The  values  may  also  be  pre-processed  to  make  the  input  signal  analytic  which  reduces  the 
number  of  cross  terms  in  the  final  distribution. 


no 


B.l  Function  myspct 


As  mentioned  in  the  text,  a  spectrogram  is  also  known  as  a  short  term  Fourier  transform 
and  is  generated  by  centering  a  window  function  at  the  time  of  interest  and  taking  a  regular 
Fourier  transform  of  the  windowed  data. 


function  spct=myspct  (a,hl)  ; 

7.  spct=niyspct (a,hl)  Calculate  the  two-sided  spectrogram 
7,  for  the  given  vector  (real  or  complex)  . 

7  If  the  transform’s  half  length  (to  be 

7  consistant  with  mybtfd  and  mywtfd)  it 

7  defaults  to  half  the  vector  length. 

7  Values  returned  are  arranged  so  the  rows 

7  correspond  to  frequencies  from  -hi  to  hl-1,  and 

7  columns  correspond  to  time  bins  (0  to  N-1) . 


7  cjm  -  31jan96 

[n,m]=size  (a)  ; 
if  n*m  >  max(Cn  m]) 

spct=’ myspct  error  -  only  works 
else 

if  n==l;  n=m;  m=l;  a=a.’;end 
tl=n;  if  nargin'"=l;  tl=hl*2;end 


7  n=# samples 
7  only  true  for  error 
with  a  vector’ 

7  make  input  a  column  vector 
7  and  set  transform  length 


tmp= [zeros (tl/2, 1) ;  a;  zeros(tl/2,l)]  ; 
spct=zeros(tl,n) ; 
for  i=l:n, 

spct (: ,i)=fftshift(fft (tmp(i:i+tl-l))) ; 

end 

spct=abs(spct) ."2; 

end 
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B.2  Function  myfwig 


This  routine  generates  the  Wigner  distribution  for  the  input  vector.  The  frequency 
resolution  depends  on  the  maximum  time  lag  applied  when  calculating  the  local  autocor¬ 
relation.  This  routine  works  with  the  alias-free  formulation  with  calculations  performed  in 
the  autocorrelation  domain  [21], 


function  wtfd=myfwig(a,l) ; 

7,  wtfd=inyfwig(a,l)  Calculate  the  two-sided  wigner  time-frequency 

7  distribution  for  the  given  signal  vector. 

7,  If  the  maximum  lag  (1)  is  not  specified,  it 

7,  defaults  to  the  vector  length. 

7,  Values  returned  are  arranged  so  the  columns 

7  correspond  to  frequencies  from  -1  to  1-1,  and 

7  rows  correspond  to  time  bins  (0  to  N-1) . 


7  this  routine  combines  code  from 
7  -  also  parallels  mybtfd.m 

7  cjm  -  27aug95  -  modified  2apr96 

[m,n]=size(a) ;  7 

if  n*m  >  max([n  m])  7 

tfd=’mybtfd  error  -  only  works 
else 

if  n==l;  n=m;  m=l;  a=a.’;end  7 
ml=n;  if  nargin"'=l;  ml=l;end  7 


my  lac. m  and  mywig.m 

to  save  memory 

n=#samples,  m  will  be  used  for  lag 
only  true  for  error 
with  a  vector^ 

make  input  a  row  vector 
and  set  maximum  lag  value 


7  first  find  the  local  autocorrelation 

maxd=min(  Cn,ml]  )  ;  7  limit  loop  to  range  of  delays 

ac=conj (a) ;  7  conjugated  a 

a=[a  zeros(l,maxd)] ;  7  padded  a 

lacm=zeros(maxd,n) ; 
for  row=l:maxd; 

lacmCrow, : )=a(row:row-i+n) .*ac; 
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end 


•/,  next  apply  the  wigner  kernel 
7.  we  know  [ml,n]  =size (lacm)  ; 
for  r=l:maxd/2; 
row=r*2; 

lacmCrow-l,  :)  =  [zeros (l,r*~l)  lacin(row~l ,  1  :n-r+l)]  ; 
lacmCrow, :)=0.5*(  [zeros(l,r-l)  lacmCrow, 1 :n-r+l)]  . . . 

+ [zeros (l,r)  lacmCrow, 1 in-r)] ) ; 

end; 

7,  finally,  convert  from  autocorrelation  domain  to  time -frequency  domain 
wtfd=real(fft (  [lacm;  zeros(l,n) ;  conj (f lipud(lacm(2:ml, :)))])); 

7i  (put  negative  frequencies  first) 

wtfd=[wtfd(ml+l:2*ml, :);  wtfd(l:ml, :)] ; 

end 
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B.3  Function  myfbin 


This  routine  generates  the  binomial  distribution  for  the  input  vector.  The  frequency 
resolution  depends  on  the  maximum  time  lag  applied  when  calculating  the  local  autocor¬ 
relation.  This  routine  works  with  the  alias-free  formulation  with  calculations  performed  in 
the  autocorrelation  domain  [21]. 

function  btfd=myfbin(a,l) ; 

y,  btfd=inyfbin(a,l)  Calculate  the  two-sided  binomial  time -frequency 

y  distribution  for  the  given  signal  vector, 

y.  If  the  maximum  lag  (1)  is  not  specified,  it 

y,  defaults  to  the  vector  length, 

y.  Values  returned  are  arranged  so  the  columns 

y,  correspond  to  frequencies  from  -1  to  1-1,  and 

y,  rows  correspond  to  time  bins  (0  to  N-1)  . 

y  this  routine  combines  code  from  mylac.m  and  mybin.m 
y,  -  also  parallels  mywtfd.m 

y  cjm  -  23aug95  updated  to  reduce  memory  use  2apr96 

[m,n]=size(a)  ;  */*  n=#samples,  m  will  be  used  for  lag 

if  n*m  >  max([n  m]  )  y  only  true  for  error 

tfd=’mybtfd  error  -  only  works  with  a  vector’ 
else 

if  n“l;  n=m;  m=l;  a=a.’;end  y  make  input  a  row  vector 
ml=n;  if  nargin"=l;  ml=l;end  y  and  set  maximum  lag  value 

y  first  find  the  local  autocorrelation 

maxd=min(  [n,ml]  )  ;  */,  limit  loop  to  range  of  delays 

ac=conj  (a)  ;  */,  conjugated  a 

a=[a  zeros (1, maxd)]  ;  */♦  padded  a 

lacm=zeros(maxd,n)  ; 
for  row=l:maxd; 

lacmCrow, : )=a(row:row-l+n) .*ac; 
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end 


um 

%  At  this  point,  lacm  holds  half  the  autocorrelation  matrix,  with  all 
7,  values  slid  to  the  left.  This  means  the  second  row  should  be  shifted 
7,  one  HALF  column  to  the  right,  the  second  row  should  shift  one  full 
7®  column  to  the  right,  etc.  This  shifting  will  be  accounted  for  when 
7o  the  kernel  is  applied  below. 

7.7®7.7.7. 


7®  next  apply  the  binomial  kernel 
7#  we  know  Cml,n]=size(lacm)  ; 
krn=  [ . 5  . 5] ; 
for  row=2:maxd; 

lacm (row, :) =conv( lacm (row, l:n-row+l) ,krn) ; 
krn=conv (krn , [ . 5  . 5]  ) ; 
end; 


7®7®7.7®7® 

%  The  matrix  lacm  now  contains  values  ready  for  transforming  via  the 
7®  FFT  into  the  time/frequency  domain.  The  values  should  correspond 
7®  to  the  upper  half  of  the  matrix  generated  via  use  of  dos_auto  and 
%  bin.m  in  the  FTP  distributed  files. 

7®7.7.7®7® 


finally,  convert  from  autocorrelation  domain  to  time -frequency  domain 
btfd=real(fft ( [lacm;  zeros (l,n) ;  conj (f lipud(lacm(2 :ml , :)))])); 

7®  (put  negative  frequencies  first) 

btf d= [btf d(ml+l : 2*ml , : ) ;  btf d( 1 :ml , : )] ; 

end 
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B.4  Function  renyi 


This  function  finds  the  Renyi  entropy  for  a  two-dimensional  array  [44].  These  function 
served  as  the  fitness  function  for  the  genetic  algorithm  kernel. 

function  h=renyi(cs,a,dt ,df ) 

7o  h=renyi(cs,a,dt ,df )  Calculate  the  Renyi  entropy  for  the  time- 
y.  frequency  distribution  cs  (matrix)  using 
y.  default  order  (a)  of  3,  and  default  delta 
X  values  (dt,  df)  of  1  (for  0  bits) 

y,  cjm  -  10oct95 

if  nargin<4;  df=l;  end  7  handle  default  values 
if  nargin<3;  dt=l;  end 
if  nargin<2;  a=3;  end 

h  =  log(  sum(sum(  (cs/sum(sum(cs) ) ) . "a  ))  *  dt  *  df)  /  (log(2)  *  (l~a)); 
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B.5  Function  genalg 


This  genetic  algorithm  tries  to  find  a  Cohen’s  class  distribution  that  satisfies  most  of 
the  desirable  properties  outlined  in  chapter  2  and  maximizes  the  Renyi  entropy  of  the  time- 
frequency  distribution.  The  routine  uses  integer  coding  of  the  genes  and  keeps  the  top  half 
of  each  generation  to  serve  as  the  parents  for  the  next  generation  [20,  34]. 

f unc t  ion  [mg , mx]  =genalg  ( si g ,  ml , nw , bw ,  pop , mr ,  mi ) 

%  [mg,mx]=genalg(sig,ml,nw,bw,pop,mr) 

•/,  genetic  algorithm  to  find  the  ’best’  kernel  for  a  Cohen’s 

y,  class  TFD.  The  signal  (sig)  goes  through  a  transform  with 

7,  max  lag  ml  and  a  kernel  using  a  cosine  series  with  nw  odd 

•/o  harmonics,  with  each  parameter  weight  coded  in  bw  bits. 

•/i  The  gene  pool  is  pop  members  long.  The  algorithm  will 

7  iterate  until  the  maximum  performer  is  constant  for  mr 

y,  cycles  or  mi  cycles  are  completed. 

y«  Returned  values  are  the  ’best’  gene  per  iteration  and  the  scores, 
y,  cjm  *“  30jun96 


iter=0; 
rpts=l; 
gl  =  nw*bw; 

mx  =  []  ;  mg  =  []  ;  lb=  []  ; 
pool  =  f loor(rand(pop,gl)+.5) ; 
score  =  zeros (pop, 1) ; 
for  i=pop/2+l :pop 

score (i) =real (renyi (myf gen (sig 

end 

while (  (iter<mi)  &  (rpts<mr)  ), 
for  i=pop/2:"l:l 
score (i)=0; 
for  pt=i+l:pop 

if  prod(pool(i, :)==pool(pt, 
end 


y,  iteration  count 
7p  number  of  repeats 
7,  gene  length 

y,  zero  out  the  ’best’  arrays 
7.  randomize  the  gene  pool 
7o  clear  out  scores; 

y,  after  first  time,  we’ll  know  these 
ml, pool (i, :) ,nw,bw))) ; 


))==!,  score (i)=score (pt) ; end 
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if  score (i)==0, 

score(i)=real(renyi(myfgen(sig,ml,pool(i, :) ,nw,bw))) ; 

end 

end 

[score, i]=sort (score)  ;  7.  arrange  by  lowest  to  highest  entropy 

pool=pool(i, : ) ; 

inx=[mx;  score  (pop)];  7#  keep  info  on  best 

nig=[mg;  pool  (pop, :)]  ; 

rpts=rpts+l; 

if  prod(lb==pool(pop, :))==0,  rpts=l;  end 

lb=pool(pop, :) ; 

iter=iter+l; 

for  i=l:2:pop/2  %  propagate  next  generation 

cp=floor(rand(l)*(gl-2)+2) ;  7.  pick  crossover  point 

pool(i  ,:)  =  [pool(pop+l-i,l:cp-l)  pool(pop/2+i ,cp:gl)]  ; 
pool(i+l,:)  =  [pool(pop/2+i,l:cp’-l)  pool(pop+l-i,cp:gl)]  ; 
end 

cp=f  loor(rand(l)*gl+l)  ;  7p  pick  random  mutation  point 

i=f loor(rand(l)*pop/2+l)  ;  7o  limit  to  new  offspring 

pool(i,cp)  =  abs(pool(i,cp)-l) ;  7o  flip  the  bit 
[iter  rpts  mx^]  7o  optional  feedback  to  user 

end 
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B.6  Function  myfgen 


This  function  generates  a  time-frequency  distribution  using  the  kernel  as  represented 
by  the  genetic  algorithm.  The  returned  matrix  exactly  matches  the  layout  generated  by 
myspct.m,  myfwig.m,  and  myfbin.m. 

function  tfd=myfgen(a,l,code,nw,bw) ; 

7«  tfd=myfgen(a,l,code,nw,bw) 

y.  Calculate  the  time -frequency  distribution  using 
7#  a  weighted  cosine  kernel.  The  n  weights  are  coded 
7.  values  between  zero  and  1  using  bw  bits.  The 
7.  weights  apply  to  an  expansion  of  odd  cosine 
70  harmonics. 

7o  Values  returned  are  arranged  so  the  columns 
7«  correspond  to  frequencies  from  -1  to  1-1,  and 
7  rows  correspond  to  time  bins  (0  to  N-1)  . 

7o  this  routine  combines  code  from  mylac.m  and  mybin.m  -  also  parallels 
70  mywtfd.m  -  derived  from  myfbin.m 

70  cjm  -  23aug95  updated  to  reduce  memory  use  2apr96 

70  cjm  -  29jun96  recoded  for  genetic  decoding  and  weighting 

[m,n]=size(a)  ;  7o  n=#samples,  m  will  be  used  later  for  lag 
if  n*m  >  max([n  m]  )  7,  only  true  for  error 

tfd=’mybtfd  error  -  only  works  with  a  vector’ 
else 

if  n==l;  n=m;  m=l;  a=a.’;end  7o  make  input  a  row  vector 
ml=l; 

7.  first  find  the  local  autocorrelation 

maxd=min(  [n,ml] ) ;  7#  limit  loop  to  interested  range  of  delays 
ac=conj  (a)  ;  7.  conjugated  a 
a=[a  zeros(l,maxd)]  ;  7.  padded  a 
lacm=zeros(maxd,n)  ; 
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for  row=l :maxd; 

lacmCrow, : )=a(row:row-l+n) .*ac; 

end 


um 

*/*  At  this  point,  lacm  holds  half  the  autocorrelation  matrix,  with  all 
y,  values  slid  to  the  left.  This  means  the  second  row  should  be  shifted 
y,  one  HALF  column  to  the  right,  the  second  row  should  shift  one  full 
y,  column  to  the  right,  etc.  This  shifting  will  be  accounted  for  when 
y,  the  kernel  is  applied  below. 

vmi 


y»  new  code  for  genetic  coding 

y. 

bwts  =  [-1  2."(-l:~l:-bw“l)]  ;  */.  numeric  weights  for  each  code  bit 

mf  =  l:2:nw*2;  */.  multiplier  for  odd  harmonics 

wts  =  zeros(nw,l);  %  find  weight  values 

whos; pause 

for  i=i:nw 

wts(i)=code((i-l)*bw+l :i*bw)*bwts;  %  code  must  be  row  vector 


end 


for  row=2:maxd  */#  now  apply  kernel  to  autocorrelation 
krn=cos(  (  (l:row) ’*pi/(row+l)~pi/2  )  *  mf  )  *  wts; 
krn=krn/sum(krn)  ;  make  sure  row  sums  to  one 
lacmCrow,:)  =  conv (lacm (row, 1 :n-row+l) ,krn) ; 

end 

y. 


y,  end  of  new  code 


y.y.y;/pyoyay.yoy.y.y.y.y.y.y.y.y. 


yxm 

y.  The  matrix  lacm  now  contains  values  ready  for  transforming  via  the 
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7,  FFT  into  the  time/frequency  domain.  The  values  should  correspond 
7.  to  the  upper  half  of  the  matrix  generated  via  use  of  dos_auto  and 
%  bin.m  in  the  FTP  distributed  files. 

7.7.7.7.7. 

7.  finally,  convert  from  autocorrelation  domain  to  time-frequency  domain 
tfd=real(fft([lacm;  zeros(l,n);  conj (flipud(lacm(2:ml, :)))])) ; 

7,  (put  negative  frequencies  first) 

tfd=[tfd(ml+l:2*ml, :) ;  tfd(l:ml, :)]  ; 

end 
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